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Abstract 


Two-sided  orthogonal  decompositions  (TSOD)  have  been  essential  tools  for  es¬ 
timating  the  numerical  rank  of  a  matrix  and  computing  various  important  subspaces 
including  the  range  (signal)  and  null  (noise,  error)  spaces.  They  include  partial  and 
complete  singular  value  decompositions  (SVD),  the  URV  decomposition  (URVD),  and 
the  ULV  decomposition  (ULVD). 

Given  the  TSOD  of  an  m-by-n  (m  >  n)  matrix  A,  it  is  often  desirable  to  suc¬ 
cessively  add  a  new  row  to  A  and  to  compute  the  TSOD  of  the  modified  matrix.  This 
is  called  the  updating  problem.  The  opposite  computation,  the  downdating  problem, 
deletes  the  existing  row  from  A,  and  computes  the  TSOD  of  the  modified  matrix.  These 
problems  of  updating  and  downdating  can  be  transformed  into  those  of  modifying  a 
symmetric  positive  definite  matrix  by  a  rank-one  matrix. 

In  this  thesis  we  propose  several  algorithms  for  rank-one  updates  and  downdates 

to  these  decompositions  with  strong  stability  properties  and  efficient  implementations 

2 

on  high-performance  computers.  We  seek  algorithms  which  only  require  0(n  )  opera- 
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tions  per  update  or  downdate  unlike  recomputing  the  TSOD  in  0(n  ).  We  also  desire 
highly  regular  data  movement  inherited  in  these  algorithms  in  order  to  implement  these 
efficiently  on  the  distributed-memory  MIMD  multiprocessors.  The  algorithms  are  based 
upon  “chasing”  strategies  for  updating  and  downdating  procedures  for  orthogonal  de¬ 
compositions. 

In  modifying  the  SVD  and  partial  SVD,  our  algorithms  separate  singular  values 
into  “large”  and  “small”  sets  and  then  obtain  an  updated  bidiagonal  form  with  corre¬ 
sponding  “large”  and  “small”  columns.  This  makes  more  accurate  update  or  downdate. 


The  algorithm  can  be  implemented  almost  identically  for  both  updating  and  downdat¬ 
ing  by  reducing  the  problems  to  a  2  x  2  updating/downdating  problem.  Moreover,  the 
bidiagonal  reduction  phase  is  highly  parallelizable.  A  perturbation  theory  for  modifying 
the  SVD  is  also  presented;  it  shows  that  the  computed  subspaces  associated  with  large 
and  small  singular  values  are  as  accurate  as  can  be  expected. 

An  alternative  to  performing  the  singular  value  decomposition  is  to  factor  a  matrix 

into  A  =  U  ^  ^  ^  VT  where  U  and  V  are  orthogonal  matrices  and  C  is  a  lower  triangular 
matrix  indicating  a  separation  between  two  subspaces  by  column  size.  These  subspaces 
are  denoted  by  V  =  (V  V  ),  where  the  columns  of  C  are  partitioned  conformally  into 
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C  =  (C^  Cy)  with  HC^IIjj,  <  e  .  Here  c  is  some  tolerance.  In  recent  years,  this  has  been 
called  the  ULVD.  A  downdating  algorithm  is  proposed  which  preserves  the  structure 
in  the  downdated  matrix  C  to  the  extent  possible.  Strong  stability  results  have  been 
proven  for  these  algorithms  based  on  a  new  perturbation  theory.  When  C  is  given  as  an 
upper  triangular  matrix,  we  have  the  URVD.  We  describe  algorithms  for  modifying  the 
URVD,  and  make  comparison  with  our  algorithms  for  modifying  the  ULVD  in  terms  of 
the  computed  subspaces. 

When  downdating  the  ULVD,  a  deflation  step  is  necessary  to  compute  its  nu¬ 
merical  rank.  We  propose  an  improved  algorithm  which  almost  always  guarantees  the 
rank-revealing  structure  of  the  decomposition  after  a  downdate  without  the  deflation 
process.  This  requires  some  condition  estimation.  Moreover,  one  can  monitor  the  con¬ 
dition  of  the  downdating  problem  by  tracking  exact  quantities  of  Frobenius  norms  of  all 
three  blocks  of  the  lower  triangular  factor  in  the  decomposition.  The  algorithm  is  also 
used  to  update  the  ULVD  with  a  slight  modification. 
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A  fully  parallel  algorithm  for  modifying  the  SVD  is  also  presented.  We  consider 
both  cyclic  and  consecutive  storage  schemes.  We  will  show  that  the  latter  scheme  out¬ 
performs  the  former  on  a  coarse-grain  distributed-memory  MIMD  multiprocessor  mainly 
due  to  high  communication  cost  required  by  the  former.  We  present  the  experimental 
results  on  the  32-node  Connection  Machine  (CM-5).  A  speed-up  of  20  and  the  efficiency 
of  CPU  utilization  60%  are  achieved  for  matrices  of  moderate  size. 

Our  algorithms  for  modifying  the  TSOD  offer  a  promising  approach  to  a  number 
of  problems  like  the  recursive  total  least  squares,  linear  regression,  the  subspace-based 
methods  for  signal  processing,  image  processing,  and  pattern  recognition.  These  prob¬ 
lems  require  a  real-time  solution  in  estimating  the  numerical  rank  of  the  data  matrices, 
and  orthonormal  basis  for  the  subspaces  associated  with  large  and  small  singular  values. 
Our  algorithms  are  capable  of  providing  those  answers  since  continual  updating  and 
downdating  are  required  by  the  underlying  physical  model. 
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If  M  is  a  partially  reduced  bidiagonal  matrix  of  the  form, 


Ik  n—k 


B1  0 


0  B2  J  n~k 


where 


ip,  V  > P2iiF<<, 


and  one  of  B ^  and  B  is  upper  bidiagonal,  and  the  other  diagonal,  (1.1)  is  called  the 
partial  singular  value  decomposition  (partial  SVD)  described  by  Van  Huffel  [118].  Here, 
we  presume  they  are  decoupled,  namely,  (&,  k  +  1)  entry  of  M  is  zero. 

If  M  is  an  upper  triangular  matrix  of  the  form, 


k  n—k 


R  S 


0  T  /  n-k 


where 


ll(sT  rT)||f.  <  t,  p ‘iij 1  >  <oi, 


and  R  and  T  are  upper  triangular,  (1.1)  is  called  the  URV  decomposition  (URVD). 
If  M  is  a  lower  triangular  matrix  of  the  form, 


k  n—k 


L  0 


F  G  I  n-k 
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where 

ll(F  G)||f  <  €,  l|i“'|f1  >  tol, 

and  L  and  G  are  lower  triangular,  (1.1)  is  called  the  ULV  decomposition  (ULVD). 

Here  k  is  the  numerical  rank  of  A  and  e  <  \Jn  —  k  *  tol.  We  use  J|  •  ||  to  denote 
the  Euclidean  norm  and  ||  *  ||^  to  denote  the  Frobenius  norm  of  a  matrix. 

The  SVD  has  been  one  of  the  most  important  tools  widely  used  in  a  number  of 
fields  of  science  and  engineering  for  decades,  mainly  because  it  offers  abundant  informa¬ 
tion  about  the  matrix  in  question.  The  SVD  has  many  benefits.  It  provides  orthonormal 
basis  for  important  subspaces  associated  with  the  matrix  including  the  range  (signal) 
and  null  (error,  noise)  spaces. 

The  URVD  and  ULVD  are  particular  cases  of  what  Lawson  and  Hanson  [70]  called 
HRK  decompositions.  Both  URVD  and  ULVD  were  introduced  by  Stewart  [101,  102],  as 
an  alternative  to  the  accurate  but  expensive  SVD.  Stewart  also  gave  methods  to  update 
these  decompositions  in  0(n  )  operations. 

In  fact,  all  of  these  decompositions  provide  valuable  information  about  the  data 
matrices.  Most  importantly  for  many  applications,  they  provide  the  orthonormal  basis 
for  the  range  and  null  spaces.  For  instance,  if  V  is  partitioned  according  to 


V 


=  (Vj  v2),  v1  e  nnxk,  v2  e  ftnx(”  k) 


(1.6) 


then  it  is  not  difficult  to  see  that  the  columns  of  V  give  the  desired  orthonormal  basis 
for  the  approximate  null  space. 
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1.1.2  Modifying  the  TSOD 


We  are  interested  in  computing  the  TSOD  of  A  when  the  TSOD  of  A  is  known, 
where  for  updating , 


A  = 


'a' 


\f  J 


(1.7) 


and  for  downdating, 


A  = 


l  T 

T 


(1.8) 


Here,  we  assume  appending  a  new  row  to  A  to  be  the  last  row  of  A,  and  deleting  the 
first  row  of  A  when  downdating.  The  downdating  problem  is  considered  more  sensitive 
than  the  updating  problem  because  small  singular  values  of  A  tend  to  diminish  after 
a  downdate,  leaving  the  matrix  near  singular,  and  thus  can  be  unstable  [99].  On  the 
other  hand,  updating  increases  all  its  singular  values.  Clearly,  refactoring  the  whole 
decomposition  without  using  the  TSOD  of  A  is  not  practical;  it  requires  0(n  )  operations 
to  compute  any  TSOD. 


1.2  Problem  Formulation 

We  transform  the  updating/downdating  problem  into  the  rank-one  modification 
of  the  symmetric  eigenvalue  problem.  Since  from  (1.7)  and  (1.8), 


-T  - 
A  A 


T  T 

A  A  +  prr 


T  T 
VM1  MV 1 


L  T 

+  prr 


V(MTM  +  pzzT)V 


T 
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where  p  >  0  for  updating  and  p  <  0  for  downdating,  and 


Thus,  the  problem  of  modifying  the  TSOD  of  A  is  equivalent  to  the  following 

T 

eigenvalue  problem:  given  a  symmetric  positive  definite  matrix  A  A  with  known  eigen- 

'T'  T  jP  T  jT 

system  A  A  =  VM  MV  ,  compute  the  eigensystem  of  M  M  +  pzz  ,  that  is,  to  find 

an  orthogonal  matrix  V  €  TZnxn  such  that 


MTM  +  pzzT  =  VMTMVT. 


(1.10) 


However,  we  do  not  form  the  explicit  product  A  A  because  of  possible  loss  of 

.  ji 

information  in  forming  A  A.  Furthermore,  the  eigendecompositions  do  not,  in  general, 
preserve  the  block  structure.  Instead,  we  compute  orthogonal  matrices  U  G 


V  G  'RnXn  such  that 


= 


for  updating, 


(1.11) 


zT?  ]  VT  =  {  ° 


for  downdating, 


(1.12) 


where  M  is  bidiagonal  for  the  (partial)  SVD,  upper  triangular  for  the  URVD,  and  lower 


triangular  for  the  ULVD. 


1,3  Importance  of  the  Problem 

Updating  and  downdating  are  important  in  signal  processing  and  statistical  appli¬ 
cations  as  new  observations  are  added,  and  the  old  observations  are  successively  deleted. 
They  can  efficiently  be  applied  to  problems  which  arise  in  a  number  of  applications: 
recursive  total  least  squares  problems  [23,  34],  linear  regression  [112,  123],  linear  predic¬ 
tion  [113],  pattern  recognition  [16],  system  identification  [60,  105],  spectral  estimation 
[17,  28,  68],  adaptive  beamforming  [76,  77,  96],  image  processing/restoration  [5,  7,  82], 
adaptive  filtering  [64],  direction  finding  [2,  73],  subspace-based  algorithms  in  signal  pro¬ 
cessing  such  as  MUSIC  (Multiple  Signal  Classification)  [94,  95]  and  ESPRIT  (Estimation 
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of  Signal  Parameters  via  Rotational  Invariance  Techniques)  [91,  92],  and  ocean  acoustic 
tomography  [110]. 

1.4  Issues  and  Concerns 

Algorithms  for  modifying  the  TSOD  must  have  at  least  the  following  features: 

Efficiency  The  algorithm  should  require  as  few  operations  as  possible,  for  example, 
0(n  ).  This  feature  would  make  it  possible  to  implement  the  algorithms  for  appli¬ 
cations  that  require  a  real-time  processing,  where  continual  updating/downdating 
of  the  decompositions  is  required. 

Stability  The  algorithm  should  produce  correct  answers  within  the  uncertainties  of 
the  given  data.  Therefore,  the  computed  solution  should  be  as  good  as  our  data 
warrants. 

Parallelism  It  should  be  easy  to  implement  the  given  algorithm  on  a  parallel  processor. 
It  is  desirable  to  develop  parallel  procedures  which  achieve  the  best  possible  load 
balance  and  minimize  the  communication  cost,  showing  high  efficiency  and  a  good 
speed-up  even  for  small  sized  problems. 

1.5  Basic  Approach  to  the  Problem 

Our  approaches  to  modifying  the  TSOD  use  ideas  from  “chasing”  algorithms 
[1,  93,  115,  125]  and  from  the  downdating  algorithm  due  to  Saunders  [46,  85].  Chasing 
algorithms  apply  a  series  of  Givens  rotations  from  both  sides  to  annihilate  all  components 
of  z,  reducing  M  to  a  desired  form.  These  algorithms  offer  highly  regular  data  movement, 
so  powerful  pipelining  strategies  can  be  used  on  parallel  computers.  A  systolic  array 
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implementation  of  a  scheme  with  similar  data  movement  patterns  to  this  one  (but  not 
similar  numerical  properties)  is  implemented  in  [117]. 

The  alternative  to  chasing  algorithms  for  modifying  the  SVD  is  that  of  finding 
the  zeroes  of  a  particular  spectral  function  [10,  25,  49,  53,  54,  67,  97].  However,  that 
approach,  as  yet,  does  not  allow  us  to  separate  the  singular  values  into  separate  blocks 
in  the  manner  discussed  in  Chapter  6. 

1.6  Main  Results 

The  following  are  the  main  results  of  this  thesis. 

•  Blockwise  procedures  for  modifying  the  TSOD  which  preserves  the  separation  be¬ 
tween  subspaces  associated  with  the  “large”  and  “small”  singular  values. 

•  An  error  analysis  of  these  procedures  demonstrating  that  the  subspaces  of  the 
modified  matrix  are  as  good  as  can  be  expected. 

•  Efficient  parallel  implementation  for  modifying  the  TSOD  that  incorporates  clever 
pipelining  strategies  using  highly  regular  data  movement  inherited  from  the  chasing 
algorithms. 

1.7  Review  of  Related  Work 

As  mentioned  in  Section  1.2,  problems  of  modifying  the  TSOD  can  be  viewed  as 
modifying  the  symmetric  positive  definite  matrix  followed  by  a  rank-one  matrix.  Gill, 
Golub,  Murray,  and  Saunders  [46]  considered  a  problem  of  modifying  the  decomposition 
of  a  matrix  following  a  rank-one  change,  where  they  showed  how  to  construct  recurrences 
for  the  product  of  Givens  rotations  in  order  to  modify  Cholesky  factor.  An  algorithm 
by  Saunders  was  also  described  for  downdating  QR  factorization,  which  has  become 
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a  backbone  of  a  number  of  downdating  algorithms  including  those  proposed  in  this 
dissertation. 

Bunch,  Nielsen  and  Sorensen  [26]  studied  a  problem  of  rank-one  modification 
of  the  symmetric  eigendecomposition  that,  in  turn,  gave  a  rise  to  their  algorithm  for 
updating  the  SVD  [25].  Their  method  was  based  on  solving  the  secular  equations. 

Dongarra  and  Sorensen  [38]  proposed  an  algorithm  that  always  computes  the 
eigenvalues  of  tridiagonal  matrices  with  high  relative  accuracy.  However,  when  eigen¬ 
values  are  clustered  together,  their  algorithms  had  difficulties  in  computing  numerically 
orthogonal  eigenvectors. 

An  improved  version  was  proposed  by  Sorensen  and  Tang  [97],  in  which  they 
incorporated  simulated  extended  precision  to  overcome  the  difficulties  of  previous  al¬ 
gorithm.  However,  using  the  simulated  extended  precision  made  the  algorithm  require 
IEEE  floating-point  arithmetic. 

With  careful  rearrangement  of  computations  in  solving  secular  equations,  Gu  and 
Eisenstat  [54]  have  succeeded  in  developing  a  backward  stable  algorithm  which  com¬ 
putes  numerically  orthogonal  eigenvectors  without  using  the  simulated  extended  preci¬ 
sion.  They  also  observed  that  by  using  the  fast  multipole  method  of  Carrier,  Greengard, 

2 

and  Roklin  [29,  30],  eigenvectors  can  be  computed  in  G(n  )  operations  as  compared 

3 

to  0(n  )  for  the  QR  algorithms  [47,  48].  They  applied  this  technique  further  to  sym¬ 
metric  tridiagonal  eigenproblems  [56],  bidiagonal  SVD  problem  [55],  and  the  problem  of 
downdating  the  SVD  [57]. 

To  this  end  several  chasing  strategies,  which  originated  from  Rutishauser  [93], 
have  been  proposed  for  updating  the  SVD  [27],  reducing  bordered  band  matrices  [1,  52, 
115],  and  their  parallel  versions  [116,  117],  and  two-way  chasing  scheme  by  Zha  [125]. 
Zha’s  algorithm  improved  conventional  one-way  chasing  procedures  by  about  a  factor  of 


2. 
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These  chasing  schemes  implemented  a  number  of  algorithms  for  modifying  the 
SVD  and  partial  SVD  [1, 14],  updating  the  URVD  [101]  and  ULVD  [102],  downdating  the 
URVD  [87]  and  ULVD  [13],  refinement  techniques  for  the  URVD  and  ULVD  [104,  100], 
and  modifying  the  ULLV  decomposition  of  two  matrices  [22,  75].  Parallel  versions  of 
some  of  these  algorithms  were  also  studied  in  [81,  103,  124]. 

1.8  Overview  of  the  Dissertation 

In  the  next  chapter,  we  review  some  fundamental  concepts  from  linear  algebra 
used  throughout  this  dissertation.  After  introducing  notations  and  basic  notions,  we 
discuss  briefly  various  types  of  orthogonal  decompositions  and  their  relations  to  the 
TSOD.  Since  we  assume  that  the  initial  TSOD  is  given  for  all  of  our  algorithms,  we 
describe  methods  for  computing  the  TSOD  as  well  as  their  numerical  properties.  We 
also  describe  the  recursive  total  least  squares  (RTLS)  problems  as  a  potential  application 
for  our  algorithms. 

Chapter  3  contains  a  detailed  description  of  basic  algorithms  frequently  used  in 
the  subsequent  chapters.  Some  of  algorithms  include  those  for  computing  and  apply¬ 
ing  a  Givens  rotation,  various  chasing  algorithms,  and  the  LINPACK  [37]  downdating 
algorithm.  We  also  give  special  treatment  for  2-by-2  updating/downdating  procedures. 

Chapter  4  discusses  methods  for  modifying  the  ULVD.  First,  we  present  a  detailed 
description  of  the  ULVD  downdating  algorithm.  An  error  analysis  for  this  algorithm  is 
also  given  to  verify  that  the  accuracy  of  the  computed  subspaces  for  large  and  small 
singular  values  is  assessed.  Finally,  we  give  numerical  tests  of  our  algorithm  in  the 
context  of  the  RTLS. 

An  improved  algorithm  for  downdating  the  ULVD  is  proposed  in  Chapter  5. 
When  downdating  the  ULVD  of  a  matrix,  a  deflation  step  is  necessary  to  compute  its 
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numerical  rank.  We  propose  an  efficient  algorithm  that  almost  always  guarantees  the 
rank-revealing  structure  of  the  decomposition  after  a  downdate  without  the  deflation 
process.  This  always  requires  some  condition  estimation.  Moreover,  we  show  how  to 
track  exact  quantities  of  Frobenius  norms  of  all  three  blocks  of  the  lower  triangular  factor 
in  the  decomposition  in  order  to  monitor  the  condition  of  the  downdating  problem.  The 
algorithm  can  also  be  used  to  update  the  ULVD  with  a  slight  modification. 

In  Chapter  6,  methods  for  modifying  the  SVD  and  partial  SVD  are  introduced. 
The  main  feature  of  these  methods  is  the  ability  to  separate  the  singular  values  into  large 
and  small  sets  and  then  obtain  an  updated  bidiagonal  form  with  corresponding  large  and 
small  columns.  A  perturbation  theory  for  updating  and  downdating  the  singular  value 
decomposition  is  also  presented. 

We  present  a  fully  parallel  algorithm  for  modifying  the  TSOD  in  Chapter  7.  Both 
cyclic  and  consecutive  storage  schemes  are  considered  in  parallel  implementation.  We 
show  that  the  latter  scheme  outperforms  the  former  on  a  coarse-grain  distributed-memory 
MIMD  multiprocessor.  We  give  the  experimental  results  on  the  32-node  Connection 
Machine  (CM-5). 

Finally,  we  give  our  conclusion  and  propose  future  work  in  Chapter  8. 
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Chapter  2 

Background 

In  this  chapter  we  present  briefly  important  concepts  from  linear  algebra  fre¬ 
quently  used  throughout  the  dissertation.  For  detailed  description  of  each  subject,  refer 
to  [50,  58,  89,  98,  106,  121,  122]. 

2.1  Notation  and  Basic  Notions 

We  use  the  following  notations  throughout  this  dissertation. 

71  Set  of  real  numbers  denoted  by  lower  case  Greek  or  lower  case 

italic  if  there  is  no  confusion 

7 Zn  Set  of  real  n— vectors  denoted  by  lower  case  italic 

7£mxn  Set  of  real  m-by-n  matrices  denoted  by  upper  case  ITALIC  or 

upper  case  Greek  letters 

I  n-by-n  identity  matrix,  that  is,  /  =  (e  ),  where  e,  = 

k—l  n— k 


O  n-by-n  zero  matrix 

n 

T 

A  Transpose  of  A 

H  T 

A  Complex  conjugate  of  A 

A~^  Inverse  of  A  £  7ZnXn,  that  is,  A  ^  A  =  AA  = 

span{a1 , . . . ,  an)  {E”=1 

range(.4)  { y  E  :  y  =  Ax,  x  6  7ln}  =  span{ctj, . . . ,  a^} 

null(i4)  {a:  E  7lm  :  Ax  =  0} 
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rank(A)  dim(range(A)) 

A(A)  Set  of  eigenvalues  of  A,  that  is,  {A  6  1Z  :  Ax  =  Xx,  0  /  x  £  TZn} 
<x(A)  Set  of  singular  values  of  A 

a. {A)  t-th  singular  value  of  A  in  nondecreasing  order,  that  is,  <?.(A)  = 

2  2 


a  (A)  Largest  singular  value  of  A 
max 

<x  .  (A)  Smallest  singular  value  of  A 
min'  ' 

k(A)  condition  number  of  A,  n(A)  =  a  ( A)la  .  (A) 
v  '  v  '  max'  "  mm'  ' 

£>(■)  g(n)  =  0(f(n))  if  there  exist  constants  c  and  N,  such  that,  for 

all  n  >  N,  we  have  g{n)  <  cf(n)  [3] 

flop  A  floating  point  operation,  that  is,  the  amount  of  work  associated 

with  an  addition,  a  multiplication,  or  a  square  root 
fi  Machine  unit 

sign(ai)  x/ |x|  if  x  ^  0;  1  if  x  =  0 

Definition  2.1  (Vector  Norms).  Letx  e  Hn-  Then 


n 

f  n  A 

11*11,  =  £  K-l. 

llzll  =  max  lx. I, 

n  nOO  i  1  2  7 

h2=  x>n 

2=1 

t 

\i= 1  / 

Definition  2.2  (Matrix  Norms).  Let  A  e  TZmxn.  Then 

\\Ax\\  l  m  n  \1 

■'VssTif’ 

where  p  —  1,2,  oo. 
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We  use  ||  •  ||  to  denote  Euclidean  norm  ||  •  ||  ,  and 


to  denote  Frobenius  norm.  It 


can  be  easily  shown  that  a  (.4)  =  ||>4||. 

max 

Next,  we  devise  a  notion  of  distance  between  subspaces. 


Definition  2.3  ([50,  p.77]).  Let  W  =  (W^  W2)  and  Z  =  ( Z2)  be  orthogonal 
matrices  where  W  ,Z  £  TZnX^  and  W,Z  £  77nx^n  *).  If  S  =  range(W )  and 

1  1  Z  Z  1  1 

S2  =  range(Z1),  then  dist(SrS2)  =  ||wj ^ ||  =  sjl  -  °min(wf’zi). 

Thus,  if  sin(0)  =  dist(5j,  52)  for  some  6,  then  6  is  the  largest  angle  between  the 
two  subspaces. 

771 X  71 

Definition  2.4.  A  matrix  A  e  77  is 


diagonal 

if  a. .  =  0, 

i] 

i  #  j; 

tridiagonal 

o 

<S»  . 

II 

© 

1  *  -  il  >  i; 

upper  bidiagonal 

if  a. .  =0, 

tj 

i  >  j  or  j  >  i  +  1 

lower  bidiagonal 

if  a..  =  0, 
*7 

i  <  j  or  j  <  i  —  1 

upper  triangular 

sa 

N  . 

II 

CD 

*  >  j; 

lower  triangular 

if  a. .  =0, 

i  <  j; 

upper  Hessenberg 

if  a  =  0, 

tj 

i>  j  +  1; 

lower  Hessenberg 

if  a. .  =  0, 
*7 

i  <  j  -  1. 

Definition  2.5.  A  matrix  A  £  TZnxn  is 
symmetric  if  A?  =  A; 

positive  definite  if  x' AT x  >0,  0  ^  x  £  77”; 

T  T  n 

positive  semi-definite  if  x  A  x  >  0,  0  ^  x  £  77  ; 

T  T 

orthogonal  if  A  A  =  AA  =  /^; 

permutation  if  A  =  (e  ,  ...,e  )  where  (s  , . . .,  s  )  is 

S,  S  1  71 

1  n 

a  permutation  of  (1, . . . ,  n). 
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||  •  ||  and  ||  •  H^,  are  orthogonally  invariant  norms,  that  is,  for  any  x  G  Hn,  A  G 
7ZmXn,  and  orthogonal  matrices,  Q  G  'JZmxrn  and  Z  G  TZnXn,  we  have 

\m\  =  \\x\\,  \\QAZ\\  =  \\A\\,  \\QAZ\\f  =  \\A\\f. 

771 V  vy 

In  theory,  if  rank(A)  =  k  for  A  G  11  ,  then  we  have 


^>^2>  -->^>0,  °k+1  =  --  =  \  =  o 


where  a .  are  singular  values  of  A  defined  in  (1.2).  However,  in  practice,  a,  is  not 

%  iC~y  1 

exactly  equal  to  zero,  but  =  ^(^)>  where  M  is  the  machine  unit.  Therefore,  we 

define  the  numerical  €  —  rank  to  be  the  largest  integer  k  such  that 


a 


k 


>  € 


that  is, 


meaning  that  there  exists  an  obvious  gap  between  the  singular  values.  A  number  of 
signal  identification  problems  assume  a  significant  gap  in  the  singular  value  spectrum  of 
the  data  matrix.  Thus,  approximate  range  and  null  spaces  associated  with  the  large  and 
small  singular  values  can  be  defined  accordingly. 


2.2  Stability  of  Algorithms 


Let  T(x)  be  a  function  of  the  input  data  x. 
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DEFINITION  2.6.  An  algorithm  for  computing  T{x)  is  backward  stable  if  the  computed 
solution  T(x)  is  the  exact  solution  of  a  slightly  perturbed  problem  with  data  x. 

This  definition  is  similar  to  that  of  Bunch  [24].  Backward  stable  algorithms  are 
very  satisfactory  although  all  of  the  algorithms  proposed  in  this  thesis  for  modifying  the 
TSOD  are  not  backward  stable.  They  are  mixed  stable  as  defined  in  the  following  sense: 

Definition  2.7.  An  algorithm  for  computing  T(x)  is  mixed  stable  if  the  computed 
solution  T(x)  is  close  to  the  solution  of  a  slightly  perturbed  problem  with  data  x. 

Definition  2.7  is  used  in  the  context  of  modifying  orthogonal  decompositions  [21, 
83,  99],  and  downdating  least  squares  solutions  [20].  We  will  also  use  this  definition  of 
stability  when  we  analyze  our  algorithms  in  the  subsequent  chapters. 

2.3  Model  of  Computation 

The  machine  unit  p  is  the  smallest  number  which  satisfies 

\fl(aopb)  —  (aop6)|  <  p  |aop6|  (2-1) 

where  op  is  one  of  the  four  arithmetic  operations,  ,  x ,  -r,  and  fl(a  op  b )  is  the  floating¬ 
point  representation  of  the  exact  result  a  op  b.  We  take  the  usual  model  of  arithmetic, 
assuming  underflow/overflow  does  not  occur, 

fl(a  op  b )  =  (a  op  6)(1  +  0)  l£l  <  A1-  (2-2) 

Most  computers  take  this  model  except  for  those  without  guard  digits  such  as  Cray  for 
which 

fi(a  oP  b)  =  0(i  +  op  6(1  +  y,  iy ,  iy  <  p 
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With  this  model  we  may  obtain  unsatisfactory  results  in  addition  and  subtraction. 

We  also  require  for  our  model 

fl(y/x)  =  y/x(  1  +  0,  |£|</*,  x  >  0.  (2.3) 

2.4  Orthogonal  Factorizations 

Orthogonal  decompositions  play  an  important  role  in  a  number  of  problems  in 
matrix  computations  such  as  least  squares  and  eigenvalue  problems.  In  this  section 
we  review  two  special  orthogonal  transformations:  Householder  reflections  and  Givens 
rotations,  and  a  family  of  orthogonal  decompositions  that  can  be  computed  by  a  series 
of  application  of  these  transformations. 

2,4.1  Householder  Transformation 

A  Householder  transformation  (reflection)  of  order  n  takes  the  form 

H  =  —  2vv^ /v^v  (2.4) 

where  v  is  called  a  Householder  vector.  It  is  easy  to  verify  that  H  is  symmetric  and 
orthogonal.  The  Householder  transformation  can  be  used  to  annihilate  a  number  of 
components  of  a  vector.  For  example,  let  x  €  7£n.  If  we  choose  v  such  that 

v.  =  0,  i  =  1, ....  A:  —  1 

i 


v.  =  x . 

X  % 


i  =  k  +  1, . .  .,n 
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where 


>  =  - signal  (E"=t »?)7, 


then  we  obtain 


T 

Hx  =  (x^,  •  •  •  >  i? 0, . . . ,  0)  . 


Note  that  we  always  choose  s  so  that  s  and  x ^  can  have  opposite  sign  to  prevent  the 
loss  of  precision  in  the  computation. 


2-4.2  Givens  Transformation 

When  one  wishes  to  zero  elements  more  selectively,  Givens  transformations  (ro¬ 
tations)  do  the  best  job.  A  Givens  rotation  takes  the  form 


J(i,j,9)  = 


0  •••  0 


—  5 


0  0 


where  c  =  cos(0)  and  s  =  sin(0)  for  some  6 . 


71 

Given  a  vector  x  6  1Z  ,  if  we  choose 


x .  -x, 

s  =  — — 


’  p=\fi 


2.2 
+  x., 


c  — 
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then  we  obtain 

J(i,mT*  =  (xl . . ,(-i'°'1i+r-,,>r 

Algorithm  3.2  described  in  the  next  chapter  computes  c  and  s  without  causing  overflow 
and  underflow  in  computing  p. 

2.4.3  Rank  Revealing  QR  Factorization 

The  QR  factorization  of  an  m-by-n  matrix  A  is  given  by 


A  =  Q 


R  » 


0  /  m—n 


Q  =  («,  Q2),  Q,  e  nmxn,  q2  e 


where  Q  G  TZmxrn  \s  orthogonal  and  R  G  lZnXn  is  upper  triangular. 

If  A  has  full  column  rank,  the  columns  of  Q  spans  range(A).  However,  when  A 
is  near  rank- deficient,  and  computing  orthonormal  bases  for  range(A)  and  null(A)  is  of 
interest,  the  factorization  of  the  form  (2.5)  is  obviously  inadequate. 

It  can  be  shown  that  with  a  careful  rearrangement  of  columns  of  A  with  some 
pivoting,  one  can  produce  another  QR  factorization  which  reveals  the  numerical  rank  of 
A.  A  rank-revealing  QR  (RRQR)  factorization  of  A  G  7ZmXn  is  any  decomposition 


k  n—k 


Rn  Rn  k 

ATl  =  Q  0  J?22  rt-k 


0 


0 


m—n 
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Q  =  (Q1  Q2),  Q1  € nmxn,  Q2en 


mx(m-n) 


where  and  R ^  are  upper  triangular,  is  well- conditioned,  ||J2“  ||~  «  cr^(A), 
II ^22 II  ~  ^jb+i^^  anc*  ^  *s  a  Permutation  matrix.  Here,  k  is  the  numerical  rank  of  A . 
Note  that  ||iZ  ||  is  not  necessarily  small  compared  to  ||/2  ||.  Thus,  the  factorization 
does  not  produce  explicitly  the  approximate  null  space  of  A. 

Chan  [31],  however,  shows  that  RRQR  factorization  may  produce  the  approximate 
null  space  for  A  with  high  accuracy  when  there  is  a  large  gap  in  the  singular  value 
spectrum  of  A  although  that  may  not  be  a  practical  assumption  in  some  applications. 
His  algorithm  incorporates  several  techniques  such  as  an  efficient  condition  estimation 
and  column  pivoting.  His  algorithm  also  guarantees  an  RRQR  factorization  for  a  high- 
rank  problem. 

Using  similar  ideas,  Foster  [44]  independently  developed  a  stable  algorithm  for 
determining  the  numerical  rank  of  a  matrix  without  requiring  column  interchanges.  Hong 
and  Pan  [62]  recently  showed  that  the  permutation  n  always  exists,  and  gave  a  method 
for  constructing  such  permutation.  However,  because  of  the  high  cost  required  by  the 
procedure,  it  has  more  theoretical  than  practical  value. 


2.4.4  Complete  Orthogonal  Decompositions 


It  turns  out  that  with  an  appropriate  choice  of  a  general  orthogonal  matrix  Z  E 
7ZnXn  (considering  permutation  matrices  as  a  special  class  of  orthogonal  matrices),  one 
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can  reduce  R  in  (2.6)  even  further  so  that  ||ii  ||  becomes  small  as  well: 


k  n—k 


Rn  Ru 


AZ  —  Q  0  R  n  k  ) 

o  0  /  m~n 


z  =  (z1  z2),  z1  e  Knxk,  z2  e  Knx(n  k) 


where  ||i2111||  1  ~  cr, .(A),  and 


"h  iW- 


Then,  it  is  easy  to  see  that  columns  of  Z2  provides  the  orthonormal  basis  for  the  ap¬ 
proximate  null  space. 


2.5  Computing  the  TSOD 


2.5.1  Computing  the  Singular  Value  Decomposition 

A  standard  way  of  computing  the  SVD  of  A  €  7 zmxn  involves  two  steps:  bidiag- 
onalization  and  computation  of  the  SVD  of  resulting  bidiagonal  matrix.  The  first  step 

TTl  X  771 

requires  to  find  products  of  Householder  transformations,  U  =  U  G  7£ 

and  V  =  V. ,  •  •  •  V  ,  6  1Znxn  such  that 
1  n— 1 


U 1  AV  = 


B  n 


0  /  m—n 
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where  B  is  upper  bidiagonal.  Then,  we  compute  orthogonal  matrices  P,Q  €  Hnxn  such 
that 

B  -  PT,QT E  =  diag(cr1,...,crn). 

Thus  we  obtain 


where  X  =  U  diag(P,  I m_n)  and  Y  =  VQ. 

The  whole  process  requires  0(m  n)  flops.  The  second  phase  of  computing  the 

SVD  of  bidiagonal  matrix  can  be  done  by  QR-iteration  [47,  48]  which  is  implemented 

in  the  LINPACK  [37].  But  the  singular  values  computed  by  this  method  differ  from  the 

true  singular  values  by  at  most  p{n)  •  //  *  a  (A),  where  p(n)  is  a  moderately  growing 

m  ax 

function  of  n .  Thus,  large  singular  values  are  computed  with  high  relative  accuracy,  but 
small  ones  are  not  generally  accurate. 

Demmel  and  Kahan  [35]  developed  an  algorithm  for  computing  all  the  singular 
values  of  a  bidiagonal  matrix  to  maximal  relative  accuracy  independent  of  their  mag- 
nitudes.  Their  algorithm  implemented  in  LAPACK  [6],  is  essentially  the  QR- iteration 
incorporated  with  a  “zero-shift”,  which  is  often  faster  than  the  standard  algorithm  im¬ 
plemented  in  the  LINPACK. 

Fernando  and  Parlett  [40]  simplified  the  zero-shift  bidiagonal  algorithm  by  Dem¬ 
mel  and  Kahan  even  further  by  replacing  a  zero-shift  QR  step  with  two  steps  of  LR 
iteration  that  implement  the  qd  algorithm.  We  describe  the  qd  algorithm  in  the  next 
chapter. 

Other  methods  for  computing  all  the  singular  values  with  high  relative  accuracy  of 
a  bidiagonal  matrix  include  bisection  [15],  Rayleigh  quotient  iteration  [88].  But  they  are 
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not  competitive  in  speed  with  the  zero-shift  bidiagonal  algorithm  and  the  qd  algorithm, 
although  they  are  probably  the  most  parallelizable  algorithms  for  this  problem. 

It  is  a  well-known  fact  that  reducing  a  dense  matrix  into  bidiagonal  form  can 
introduce  large  relative  errors  in  its  singular  values.  The  Jacobi  method  for  computing 
the  SVD  of  a  dense  matrix  is  much  slower  but  more  accurate  than  any  algorithms  that 
first  bidiagonalize  the  matrix  [36].  In  this  iterative  algorithm,  a  series  of  Givens  rotations 
are  applied  to  pairs  of  rows  and  columns  to  reduce  off-diagonal  entries.  With  a  clever 
ordering  [74],  the  algorithm  can  be  implemented  in  parallel,  being  competitive  in  speed 
with  other  methods. 

2.5.2  Computing  the  ULV  and  URV  Decompositions 

The  ULVD  of  A  €  nmxn  can  be  obtained  by  computing  its  QL  factorization 


771 X  771  Tl  X  72 

where  Q  £  1Z  is  orthogonal  and  L  £71  lower  triangular,  followed  by  computing 
the  ULVD  of  L  using  the  deflation  technique  described  in  [100,  102].  First,  we  estimate 
an  approximate  left  singular  vector  u ^  of  unit  norm  of  L  which  corresponds  to  cr^(I/) 
using  some  condition  estimator.  A  survey  of  popular  condition  estimators  is  given  in 
[61].  Then,  we  compute  an  orthogonal  matrix  Q  £  TZ  such  that  u^Q  =  and  an 

-  72  X  77 

orthogonal  matrix  P  €  11  ,  such  that 


an(L)  =  ||/i||  =  \\(Jq)(QTLP)\\  =  \\eTn(QTLP)l 
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which  is  the  size  of  the  last  row  of  L.  Here,  P  is  applied  to  restore  Q  L  into  the  lower 

triangular  form.  We  repeat  this  deflation  process  until  all  the  small  rows  of  L  appear  in 

the  decomposition  yielding  the  ULVD  of  A  of  the  form  (1.5). 

Similarly,  the  URVD  of  A  can  be  obtained  by  computing  its  QR  factorization  of 

the  form  (2.5)  followed  by  computing  the  URVD  of  R  using  the  deflation  steps.  This  time 

we  estimate  an  approximate  right  singular  vector  v  of  unit  norm  of  R  which  corresponds 

—  nxn  ~T 

to  cr ^(R ),  and  then  compute  an  orthogonal  matrix  Q  £  H  such  that  Q  v^  =  and 
an  orthogonal  matrix  P  £  'RnXn ,  so  that 

<?„(*)  =  ll*«„ll  =  \\(fTRQ)(QT\)\\  =  \\(PTRQ)'J 

which  is  the  size  of  the  last  column  of  R .  Here,  P  is  applied  to  restore  RQ  into  the 
upper  triangular  form.  We  repeat  this  deflation  process  until  all  the  small  columns  of  R 
appear  in  the  decomposition  yielding  the  URVD  of  A  of  the  form  (1.4). 

2.6  Subspaces  from  the  URV  and  ULV  Decompositions 

In  this  section  we  present  error  bounds  for  accessing  the  accuracy  of  subspaces 
computed  by  the  TSOD,  particularly,  the  ULVD  and  URVD.  The  discussion  that  follows 
is  largely  abridged  from  that  of  Fierro  and  Bunch  [41,  42]. 

Let  the  SVD  of  A  £  Kmxn,  m  >  n  be 
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where 


k  n—k  m—n 


k  n—k 


U=  (P1  V2  V 3  )•  V=  (V,  V2  j, 


and  S  has  the  form  (1.2),  and  let  the  URVD  of  A  be 


A  =  U  R  VT 
A  VR  L  1 


where 


k  n—k  m—n 


k  n-k 


UR1  UR2  UR3  )’  VL  ~  (  VRl  VR2 


and  C p  has  the  form  (1.4),  and  let  the  ULVD  of  A  be 


a  =  ul  L  Vl 

l  0  , 


(2.10) 


where 


k  n—k  m—n 


Ur=  Kl  UL2  UL*l  Vr=  Kl  Vl 


k  n-k 


LI  L2 


and  C l  has  the  form  (1.5).  Here,  U,  U^,  U^,  V ,  V ^  and  are  orthogonal  matrices. 

Then  we  have  the  following  bounds  on  subspaces  computed  by  the  URVD  and 
ULVD,  which  are  associated  with  the  large  and  small  singular  values. 
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Theorem  2.1  ([42]).  Let  A  €  HmXn  have  the  SVD  of  the  form  (2.8)  and  the  URVD 
of  the  form  (2.9).  Suppose  o  .  ( R )  >  ||T||.  Then 


11*11 


imi + imi 


*  ii™^  ? 


m-imr 


(2.11) 


min 


\W^u 


121"  “  2 


11*11  imi 


(2.12) 


Theorem  2.2  ([42]).  Let  A  G  Umxn  have  the  SVD  of  the  form  (2.8)  and  the  ULVD 
of  the  form  (2.10).  Suppose  cr  .  ( L )  >  ||G||.  Then 


W* 1  ^  * 


11*11  Ill'll 


^unW-IW 


(2.13) 


- — - <  II UTU 

\\L\\  +  \\G\\-"U2U 


LI"  ~  a 


11*11 


mm 


(L)  -  ||G|| 


(2.14) 


From  these  theorems  we  immediately  realize  that  the  quality  of  subspaces  com¬ 
puted  by  the  URVD  and  ULVD  do  not  depend  on  the  gap  in  the  singular  value  spectrum, 
which  is  required  by  the  RRQR  decomposition.  We  also  observe  that  by  keeping  ||S||  and 
|| jP||  small,  we  can  obtain  highly  accurate  subspaces.  Several  ways  of  achieving  this  are 
described  in  [41,  42,  104].  It  can  be  also  argued  that  the  ULVD  computes  the  numerical 
null  spaces  more  accurately  than  the  URVD,  whereas  the  URVD  yields  a  better  estimate 
of  the  numerical  range. 
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2.7  Total  Least  Squares  Problem 

2.7.1  Problem  Formulation 

The  classical  linear  least  squares  (LS)  problem  is  given  by 

min  || Ax  -  6||,  A  e  7Zmxn,  b  e  V™  (2.15) 

xenn 

or  equivalently, 

min  ||e||,  e  G  7 Zm  (2.16) 

6+eGrange(i4)  V  7 

where  A  is  data  matrix  whose  rows  contain  measurements  from  the  model  under  con¬ 
sideration,  and  b  is  the  observation  vector.  Here,  we  presume  A  is  free  of  error  and  b 
is  subject  to  error.  However,  it  is  unrealistic  to  take  error-free  measurements  from  the 
model. 

The  total  least  squares  (TLS)  problem  assumes  that  there  are  errors  in  the  data 
matrix  A  as  well  as  the  observation  vector  b.  The  TLS  problem  has  the  formulation 


min 

6+eGrange(j4-f  E) 


ll(£  e)||F, 


(2.17) 


which  is  analogous  to  (2.16). 

The  errors-in-variables  problems  have  a  long  history  in  statistical  literature.  In 
the  field  of  numerical  analysis,  this  problem  was  first  introduced  by  Golub  and  Van 
Loan  [51]  and  then  studied  extensively  by  Van  HufFel  et  al.  [114,  118].  If  e ^  and 
( E  e  )  are  the  LS  and  TLS  corrections  to  (2.16)  and  (2.17),  respectively,  it  can 

1  1  L/O 

be  shown  that 

W(ETLS  eTLS^F  -  II eLS II 


2.7.2  Basic  Solution 


The  TSOD  gives  an  elegant  way  of  solving  the  TLS  problem.  Suppose  the  TSOD 
of  (A  b)  £  7 ^mx(n+l)  .g  gjven  by 


(A  b)  —  U 


M 

0 


k  n—k+1 


where  U  £  7^mxm  and  y  g  7j(n+1)x(n+1)  are  orthogonal,  and  M  £  7j(n+l)x(n+l)  has 
one  of  the  forms  (1.2)-(1.5).  Then  =  (t;^  j,  is  a  basis  of  the  noise  subspace 

of  (A  b),  and  the  minimum  norm  TLS  solution  is  given  by  computing  a  Householder 
transformation  H  G  72,(n+l)x(n'*'l)  such  that 


(  Y  d  \  n 

V2H  =  .  (2.18) 

v°  *)  1 

If  S  0,  the  TLS  solution  x  is  given  by 

x  =  -d/6.  (2.19) 

See  [118]  for  details.  In  particular,  Golub  and  Van  Loan  formulated  the  solution  based 
on  the  SVD.  Van  Huffel  and  Zha  [119]  also  formulated  the  solution  to  the  TLS  problems 
based  on  the  URVD  and  the  ULVD  without  the  explicit  computation  of  the  approximate 


null  space  basis  V^. 
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2.7.3  Recursive  Total  Least  Squares 

When  one  does  not  have  complete  knowledge  of  the  data  in  a  given  model,  recur¬ 
sive  procedures  make  it  possible  to  achieve  satisfactory  results.  The  algorithm  is  given, 
to  begin  with,  incomplete  knowledge  about  the  environment,  and  modifies  the  processing 
model  in  an  adaptive  fashion  as  data  are  received  sequentially. 

In  recursive  TLS  problem,  it  is  required  to  append  a  new  row  to  the  data  matrix 
A  and  an  observation  to  6,  and  the  new  information  must  be  incorporated  into  the  TLS 
solution.  It  is  also  desirable  to  delete  the  oldest  observation  from  the  existing  data. 

The  updating  and  downdating  problems  of  the  form  (1.7)-(1.8)  are  associated 
with  the  sliding  window  method  [4,  19,  32].  At  each  step  of  the  sliding  window  method 
with  the  window  size  s,  an  s  x  n  data  matrix  is  constructed  from  an  m  X  n  observation 
matrix  A  by  adding  a  new  row  to  the  data  matrix  in  the  previous  window  and  deleting 
the  oldest  row  from  it.  In  step  j,  the  row  s  +  j  of  the  observation  matrix  is  added  and 
the  row  j  is  deleted,  giving  the  data  matrix  A 


l  T 
aj-s+ 1 
T 

aj-s+ 2 


T 

a. 

J-l 


T 
a . 
J 


It  \ 

aj-s+ 2 

T 
a . 

J-s+ 3 


Aj+1 
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An  alternative  approach  is  to  nse  an  exponential  forgetting  factor  ft  (0  <  (3  <  1) 
[59].  In  this  approach  the  modified  matrix  in  (1.7) — (1.8)  is  given  by 


Here,  the  effect  of  old  observation  diminishes  exponentially  as  continuous  updating  is 
required.  However,  the  explicit  removal  of  the  observation,  as  in  the  sliding  window 
method,  makes  it  simple  to  estimate  the  rank  of  modified  matrix,  that  is,  it  remains  the 
same  or  increases  by  one  for  updating,  or  decreases  by  one  for  downdating.  Thus,  an 
indefinite  number  of  condition  estimation  steps  is  not  necessary  for  rank  detection  as 
done  in  [101,  102]. 
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Chapter  3 


Basic  Algorithms 


This  chapter  contains  detailed  description  of  basic  algorithms  frequently  used  in 
the  subsequent  chapters.  Throughout  the  dissertation  we  follow  the  convention  of  the 
Matlab  [78]  in  describing  the  algorithms.  Some  of  the  algorithms  include  those  for 
computing  and  applying  a  Givens  rotation,  various  chasing  algorithms,  and  the  LIN- 
PACK  [37]  downdating  algorithm.  We  also  give  a  special  treatment  for  2-by-2  updat¬ 
ing/downdating  procedures. 

3.1  Givens  Rotations 

This  section  contains  simple  routines  for  constructing  and  applying  Givens  rota¬ 
tions  described  in  Section  2.4.2. 

Algorithm  3.1.  This  function  computes  complex  a6solute  value  of  x  +  iy,  that  is, 
t  =  yfx^  +  y2. 

function  t  =  cabs(a:,y) 
if  \x\  >  \y\  then 


else 

r  *-  x/y;  t  <-  \y\  *  \/l  +  r2 

endif 


end 
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end 
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3.2  Chasing  Algorithms 

3.2.1  A  Chasing  Routine  for  a  Bidiagonal  Matrix 

The  following  routine,  for  a  given  vector  z  and  a  diagonal  matrix  B,  finds  orthog¬ 
onal  matrices  U  and  V  such  that 

B  =  UTBV ,  VTz  =  pev  p  -  \\z\\  (3.1) 

where  B  is  lower  bidiagonal. 

Fig.  3.1  illustrates  the  reduction  steps.  The  right  arrows  denote  rotations 
from  the  left  on  two  particular  rows,  whereas  the  downarrows  J.J  denote  rotations  from 
the  right  on  two  particular  columns.  Here  z  denotes  elements  in  the  vector  z,  b  denotes 
elements  in  B,  and  z  or  b  denotes  an  element  about  to  be  zeroed  out.  We  now  formally 
present  the  algorithm  below. 

Algorithm  3.4  (Chasing  algorithm  for  bidiagonal  reduction).  Given  a  diag¬ 
onal  matrix,  B  —  diag(7(l:  n)),  and  a  vector  to  be  reduced,  z(l:  n ),  the  following  chasing 
scheme  produces  a  lower  bidiagonal  matrix  B  such  that  B  =  bidiag(7(l:n),</>(l:n  —  1)) 
and  satisfies  (3.1). 

procedure  forchase(7,  </>,  z,n) 
formrot  (z(n  —  1),  z(n),  cn,  sn) 

e  * - sn  *  7 (n  —  1);  7(n  —  1)  <—  cn  *  ~/(n  —  1) 

< p(n  —  1)  <—  sn  *  7(n);  7(n)  «—  cn  *  7 (n) 
formrot  (7 (n),e,cn,sn) 
applyrot  (cn,  sn,  <t>(n-  l),7(n—  1),1) 
for  i  <—  n  —  2, . . .,  1 


formrot  ( z(i ),  z(i  +  1),  cn,  sn) 
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Fig.  3.L  Forward  Chasing  Procedure  for  the  Bidiagonal  Reduction 
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e< - sn  *  7(i);  ~f(i)  <—  cn  *  7(1) 

< fi(i )  <—  sn  *  7(1  +  1);  7 (i  +  1) <—  cn  *  7(1  +  1) 
d  *—  cn*  <f>(i  +  1);  <f>(i  +  1)  <—  sn  *  <f>{i  +  1) 
formrot  (7 (i),e,cn,sn) 
applyrot  (cn,  sn,  <f)(i),  7 (i  +  1),  1) 
applyrot  (cn,  sn,  d,  (f>(i  +  1),  1) 
for  j  <-  i, . . . ,  n  -  3 

formrot  (<l>{j),  d,  cn,  sn) 
applyrot  (cn,sn,7(j  +  1),  </>(j  +  1),  1) 
e  *-  sn  *  7 (j  +  2);  7 (;  +  2)  «-  cn  *  7 (j  +  2) 
formrot  (<l>(j),  d,  cn,  sn) 
applyrot  (cn,sn,7(j  +  l),<f>(j+  1),1) 
d  <—  sn  *  ^>(j  +  2);  <f>(j  +  2)  <—  cn  *  (f>(j  +  2) 
endfor 

formrot  ( <f>(n  -  2),  d,  cn,  sn) 
applyrot  (cn,  sn,7(n—  \),4>{n  —  1),1) 
e  <—  sn  *  7(n);  7(n)  <—  cn  *  7 (n) 
formrot  (7(n  -  1),  e,  cn,  sn) 
applyrot  {cn,sn,4>{n-  l),7(n),l) 
endfor 
end 


2 

This  algorithm  constructs  and  applies  n  -  n  Givens  rotations.  A  similar  proce¬ 
dure  backchase  to  that  illustrated  above  can  be  used  to  produces  orthogonal  matrices 
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U,Ve  nnxn  such  that 


B  =  UTBV,  VTz  =  pen ,  p  =  ||z|| 


(3.2) 


where  £  is  upper  bidiagonal. 

For  the  sake  of  brevity,  we  do  not  present  backchase.  It  is  computed  by  reversing 
the  two  vectors  7(1:  n)  and  (f>(l:n-  1),  performing  forchase  and  reversing  the  vectors 
back.  The  algorithm  would  just  have  the  loop  in  forchase  go  backward  instead  of 
forward. 

3.2.2  A  Chasing  Routine  for  a  Lower  Triangular  Matrix 

We  now  describe  a  simple  chasing  routine  for  a  lower  triangular  matrix.  This 
routine,  for  given  a  vector  z  and  a  lower  triangular  matrix  C ,  finds  orthogonal  matrices 
U  and  V  such  that 

C  =  UTCV ,  VTz  =  pev  p=\\z\\  (3.3) 

where  C  is  lower  triangular. 

Consider  the  4x4  case  in  Fig.  3.2,  We  state  the  procedure  lchase  formally  below. 

Algorithm  3.5.  Given  a  lower  triangular  matrix  C  and  the  updating  vector  z,  this 
procedure  performs  a  chasing  operation  on  C,  and  produces  a  lower  triangular  matrix 
C  that  satisfies  (3.3). 

procedure  lchase(c,z,n) 
for  i  <—  n  —  1, . . .,  1 

formrot  1  ),cn,sn); 

e  < - sn  *  c(i,  i);  c(i,  i)+-cn*  c(t,  i); 

applyrot(c(i  +  1:  n,i),  c(i+  1:  n,i  +  1  ),cn,  sn,  n  —  i ); 


37 


I  l  l  l 


Fig.  3.2.  Chasing  Steps  for  a  Lower  Triangular  Matrix 

formrot  ( c(i  +  1,  i  +  1),  e,  cn,sn )  ; 
applyrot(c(i  +  1, 1:  t),  c(i,  1:  i),  cn,  sn,  i ); 

endfor 

end 

A  similar  chasing  procedure  can  be  specified  for  an  upper  triangular  matrix  when 
modifying  the  URVD.  Stewart  [102]  points  out  that  if  the  matrix  C  is  from  a  rank 
revealing  decomposition  with  k  large  rows  and  n-k  small  rows,  this  algorithm  can  yield 
k  +  1  large  rows,  thus  the  rank  revealing  nature  of  C  may  be  lost. 

3.3  qd  Procedure 

It  is  easy  to  find  an  orthogonal  matrix  Q  as  a  product  of  Givens  rotations  such 

that 


B  =  QB 
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where  B  is  lower  bidiagonal  and  B  is  upper  bidiagonal.  This  is  the  same  as  one  unshifted 
Fernando-Parlett  qd  step  [40],  hence  the  name.  Fig.  3.3  illustrates  the  reduction  steps 
for  a  4  x  4  case. 


Fig.  3.3.  One  Step  of  qd  Procedure 

Algorithm  3.6.  This  procedure  produces  the  orthogonal  factorization  of  a  lower  bidi¬ 
agonal  matrix.  7(1:  n )  is  the  diagonal  on  input  and  output.  <f>(  1:  n—  1)  is  the  subdiagonal 
on  input  and  the  superdiagonal  on  output. 

procedure  qd  (7,  <f>,  n ) 
for  i  «—  1 —  1 

formrot(7(i),  (f>(i),  cn,  sn ) 

<f>(i)  <—  sn*  7  (i  +  1) 

7  (i  +  1)  <—  cn  *  7  (i  +  1) 
endfor 


end 
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3.4  The  LINPACK  Downdating  Procedure 

The  following  downdating  procedure  due  to  Saunders  [46]  is  considered  the  most 
accurate  downdating  procedure  that  does  not  require  information  from  the  first  row  of 
U  in  (1.1)  [20]  (if  we  have  that  first  row,  we  obtain  a  procedure  that  is  backward  stable 
in  the  strong  sense).  It  is  the  procedure  that  is  implemented  in  the  LINPACK  [37]. 

Algorithm  3.7  (The  LINPACK  downdating  procedure).  Given  M  G7£nXTland 

n  -  t  -  T 

zeK  ,  this  algorithm  computes  the  downdated  matrix  M,  that  is,  M  M  =  M  M  — 

T 

zz  . 

Step  1.  Solve 

M^a  =  z  (3.4) 

T  T 

If  ||a||  >  1  declare  M  M  —  zz  indefinite  and  stop.  Otherwise  go  to  Step  2. 

Step  2.  Compute  a  =  \fl  -  |)a|]2  and  Q  =  ■  •  -Qn  €  be  a  prod¬ 

uct  of  Givens  rotations,  Q .  =  J(l,t+  1 ,  ^. ) ,  i  =  1  ,...,n  such  that 


Step  3.  Compute 


We  note  that  if  M  is  upper  or  lower  triangular,  it  is  simple  to  choose  Q  as  a  product 
of  Givens  rotations  Q  ,Q^, . .  .,Qn  so  that  M  remains  upper  or  lower  triangular.  Pan 
[85]  shows  that  for  M  upper  triangular,  this  method  can  be  sped  up  by  combining  the 
forward  substitution  phase  with  the  application  of  the  Givens  rotations. 
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3.5  2x2  Updating/Downdating  Procedure 

The  following  algorithm  computes  orthogonal  matrices  G ^  —  7(1,2,^),  G  = 
7(1,3,<£  )  G  TZ^X^  for  some  <f>.,  i  =  1,2,  such  that 

£  t 


\ 

\ 

0  0 

1 

£  p 

a  b 

II 

to 

a  b 

v°  *7 

1 

\°  C  ) 

Algorithm  3.8  (2x2  Updating).  Given  scalars  a,  b,  c,  £,  p,  this  algorithm  computes 
scalars  5,  b ,  c  of  an  updated  matrix  defined  in  (3.5). 

function  [a,b,c\  =  up22  (a,b,c,£,p) 
a  =  cabs(a,  £);  a  =  a/a;  f  =  £/a 
b  =  ab  —  £p;  c  =  ap  —  £b 

end 


Similarly  we  give  an  algorithm,  which  is  based  on  Algorithm  3.7  for  downdating  a 
2x2  matrix,  that  is,  computing  orthogonal  matrices  Q  =  7(1, 3, 0j),  Q 2  =  7(1,2, 6^)  G 


such  that 


\ 

\ 

1 

£ 

P 

a 

0 

0 

0 

5 

6 

T  T 
=<«; 

a 

6 

0 

?  ) 

^7 

0 

C  ) 

(3.6) 


where  we  solve 


\ 

T 

n  \ 

,  \ 

a  b 

(3 

= 

£ 

\°  c ) 

l7  J 

\PJ 

,  a  =  \!  l-/32-72, 


(3.7) 
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where  J(i,k,6)  is  a  Givens  rotation  in  the  (z,  k)  coordinate  plane  for  some  0. 

Algorithm  3.9  (2x2  Downdating).  Given  scalars  a,  6,  c,  /3,  and  7,  this  algorithm 
computes  entries  of  downdated  matrix  defined  in  (3.6),  namely,  2,  6,  and  c. 

function  [2,6,c]  =  down22  (a,6,c,/3, 7) 

Q:  =  y/1  +  /3y/l  ~ 
a  =  y'a  +  7^/q— ^ 
p  —  ~fc/a 

c  =  ac/S;  a  =  5a;  b  =  56  —  p/3 


end 
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Chapter  4 


Modifying  the  ULV  Decomposition 


4.1  Introduction 


We  give  methods  for  modifying  the  ULVD.  Rewrite  the  ULVD  of  A  €  Tlmxn  0f 


(1.1)  as 


A  =  U 


where  U  €  'JZmxrn  and  V  6  7£nxn  are  orthogonal,  and 


k  p-k  n-p 


L  0  0 


F  G  0  P-fc 


00  0  n-p 


so  that  our  algorithm  of  ULVD  separates  out  columns  that  are  exactly  zero.  Here  C 
takes  the  position  of  M  in  (1.1). 

We  also  rewrite  z  of  (1.9)  as 


x  1  k 


z  =  VTr  = 


n—p 
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As  in  the  updating  routine  of  Stewart  [102],  the  matrices  C  and  V  can  be  produced 
using  0(n )  Givens  rotations,  thus  updating  the  factorization  in  0(n  )  operations.  Our 
approach  to  downdating  the  ULVD  (4.2)  uses  ideas  from  chasing  algorithms  [102]  and 
from  the  downdating  algorithm  due  to  Saunders  [46,  85]. 

The  following  are  the  main  results  of  this  chapter: 

1.  A  blockwise  procedure  for  downdating  the  ULVD  that  yields 

_  \ 

LOO 

C=  F  G  0  (4-4) 

^  0  0  0  y 

where 

ll(f  <5)11  <  ||(f  0)11,  ||(f  <S)||p  <  |l(f  G)||f 

and  L  and  G  are  lower  triangular,  and  the  blocks  are  conformal  with  (4.2). 

T  T  . 

2.  A  downdating  algorithm  that  works  whenever  L  L  —  xx  is  positive  definite,  the 
same  as  if  we  were  downdating  only  L.  Our  technique  maps  back  onto  the  original 
matrix  A  in  a  more  satisfactory  manner  than  the  technique  given  by  Park  and 
Elden  [87]. 

3.  An  error  analysis  of  this  procedure  showing  that  the  singular  subspaces  of  the  up¬ 
dated  matrix  are  as  good  as  can  be  expected.  We  also  give  some  new  perturbation 
results  showing  that  the  condition  of  the  downdate  is  related  only  to  the  L  block 
in  (4.4).  Thus  tracking  the  ULVD  is  a  very  stable  method  for  tracking  subspaces. 

Our  ULVD  downdating  algorithm  is  proposed  in  detail  in  Section  4.2.  Section  4.3 
gives  an  error  analysis  of  the  algorithm.  The  accuracy  of  the  computed  subspaces  for 
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large  and  small  singular  values  is  assessed.  In  Section  4.4,  we  give  numerical  tests  of  our 
algorithm  on  recursive  total  least  squares  problems. 


4.2  A  Procedure  for  Downdating  ULV  Decomposition 


4.2.1  Description  of  the  Algorithm 


We  introduce  the  following  algorithm  for  downdating  the  ULVD.  Our  procedure 

_  n~<  rp 

produces  a  downdated  matrix  C  if  for  L  and  x  defined  in  (4.2)  and  (4.3),  L  L  —  xx  is 
positive  definite. 


Algorithm  4.1  (Procedure  for  ULVD  downdating).  Given  a  lower  triangular  ma¬ 
trix  C  of  the  form  (4.2)  and  a  vector  z  of  the  form  (4.3),  this  algorithm  finds  a  lower 
triangular  matrix  C  of  the  form  (4.4)  and  orthogonal  matrices  U  and  V  satisfying  (1.12). 
Initially,  V  =  V.  The  components  in  yQ  are  ignored  (will  be  justified  by  Proposition  4.2). 
Throughout  the  description  of  this  algorithm  and  G^  denote  lower  triangular  ma¬ 
trices,  =  [F^]Tev  and  =  e^G^je^ 

Step  1.  Compute  orthogonal  matrices  U^,V^  E  such  that 


G{l)  =  U*GVv  V^y  =  pe(P-k\  p=||y||. 


Also,  compute 


P(1)  =  U^F. 


Update  U  «-  diag(Ik+vUvI '  );  V  «-  V  diag  (1^,1  ). 


45 


46 


else  <  p 


p-6p 

h 


/  a® 

frT  9n 


i<?  =  o 


f2)  T 

where  6p  —  p  —  (afifjj  +  a  h). 


endif 


Define  =  J(l,k  +  1,6)  as  the  Givens  rotation  with  cn  =  cos(0),  sn  =  sin(0) 


such  that 


9U^9U  if  9U  ^  0 

0  if  =  0. 


Update  U  <-  UU.  diag(t/  ,  /  ,). 

4  o  71 — K 


Step  5.  Find  an  orthogonal  matrix  £  7£(k+l)x(fc+l)  suc^  ^at 


(i(3)  0)  =  (Z(2)  heTx)V2 


(jp(3)  <?(<))  =  (F(2)  G(3))d iag^J^p 


Update  V  ^  Vdiag(y2,/n_it_1).  If  g^  /  0, 


- 

L 

0 

0 

F 

G 

0 

K  0 

0 

°y 

L^>  0  0 

f( 3)  G(4^  0 
0  0  0 


and  go  to  Step  7;  otherwise,  go  to  Step  6. 
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Step  6.  If  <r  /  =  0  then  / j  —  0  also  since  it  was  formed  from  ^  '  using  V  .  Thus 

J 1 1  1  1  L  Z 


k  P~k 


k  LK  ’  0 


p-k  \  J?(3) 


it  p-* 


L  0 


=  0  0 


F  G^  I  P-k-l 


where  is  a  lower  Hessenberg  matrix.  We  then  find  an  orthogonal  matrix 

V  E  7z(p~k)x(p-k)  an(j  an  orthogonal  permutation  matrix  U  E  T^P 
3  s 


such  that 


vTg^v  = 

D  u 


G  0  1  P-*-1 


0  0/  1 


Update  U  <-  U  diag^^  )\  V  «-  V  diag^,  VyIn_p).  Thus  C  has  the 


Jt  p-k- 1  n-p+1 


L  0 


C  =  f  G  0  P~k~l 

0  0  0  I  n~P+ 1 


Step  7.  Perform  a  ULVD  of  L  to  determine  its  numerical  rank.  If  we  have  deter¬ 
mined  the  rank  of  L  correctly,  the  rank  of  L  should  be  k  or  k  —  1.  Make  appropriate 
adjustments  to  F  and  G. 


We  give  an  expression  for  U  in  the  statement  of  the  algorithm,  but  our  analysis 
assumes  that  the  left  orthogonal  matrix  U  is  not  saved.  Although  it  is  not  computed  by 
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the  algorithm,  for  the  discussion  that  follows  we  need  to  define  the  vector  z  by 

(4- 

Algorithm  4.1  requires  11A;2  +  6(p  -  k)2  +  12k(p  -  k)  +  0(p)  flops  for  Steps  1-5  and 

2 

4k  +  O(k)  flops  for  Step  6.  When  V  is  modified,  additional  6 n(p  —  k)  +  6 nk  flops  will 
be  required.  The  deflation  step  requires  about  vk  +  16nk  flops,  where  v  is  the  constant 
depending  upon  the  condition  estimators  used  [61]. 

Fig.  4. 1-4.2  illustrate  the  action  of  our  algorithm  on  a  6  x  6  matrix  with  k  —  3. 
Here  Z,  g,  and  /  denote  components  of  L,  F  and  G,  and  x,  y,  and  h  are  components 
of  those  vectors.  Rightarrows  and  downarrows  denote  premultiplication  and  postmulti¬ 
plication  by  orthogonal  transformations  from  the  left  and  from  the  right,  respectively. 
Here  p,  h  and  h  have  meanings  consistent  with  those  in  the  statement  of  Algorithm  4.1. 

A  set  of  denotes  the  applications  of  and  in  Steps  3  and  4  when  downdat- 
ing  L  and  g Step  6  is  illustrated  in  Fig.  4.2.  Note  that  for  Steps  3-5,  the  components 
x  and  y  have  a  different  interpretation  from  Steps  1  and  2  above,  they  are  now  the  values 
“downdated.”  To  illustrate  their  action  properly,  U  and  U  are  given  in  reverse  order 

O  T 

rather  than  the  order  that  they  are  actually  applied. 

Remark  4.1.  This  algorithm  works  if 


Fig.  4.1.  Reduction  Steps  for  Downdating  the  ULVD 
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Fig.  4.2.  Reduction  Steps  When  G  Becomes  Singular 


T  T  T 

is  substituted  for  G  and  ( y  y Q  )  is  substituted  for  y .  That  is,  it  is  not  necessary 
to  explicitly  handle  the  zero  block,  it  can  be  made  part  of  G .  That  is  the  original 
formulation  in  [102].  However,  if  that  is  done,  whenever  y Q  ^  0,  any  zero  diagonal  will 
be  chased  to  the  g ^  position,  all  of  ( y  yQ  )  will  be  treated  as  “noise”.  If  ||G||  >  g, 
as  is  often  true  in  practice,  then  ||2/||  is  possibly  significant  whereas  yQ  can  only  result 
from  computing  errors  from  computing  (4.2)  or  multiplying  z  =  V  r.  Our  formulation 
neglects  part  of  y  only  if  the  downdate  of 


0 

G 


\ 

/ 


...  ,  T  T\T 
with  [x  y  ) 


cannot  be  done. 


In  updating,  there  is  no  similar  benefit  to  separating 


out  the  zero  block. 
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Remark  4.2.  We  note  that  (4.6)  in  Step  4  is  equivalent  to  computing 


\ 

l 

p 

to 

II 

v  ~h  ) 

OJ 

l  „(2) 


rn  +  (^>)Q 


\ 


h  + ( 8p)a 


6p  =  p-  (ag^  +  aT h) 


~T  ,  T  T  .  ~ 

since  f/  maps  the  additive  noise  vector  (<$p)(a  a  )  into  ( 6p)e y  Only  h  is  of  interest 

in  the  computation. 

C  satisfies  (4.4)  as  is  proven  quickly  in  the  following  proposition. 


Proposition  4.1.  For  the  matrix  C  resulting  from  Algorithm  4-1  we  have  (4-4)- 

Proof.  This  proof  is  a  straightforward  consequence  of  facts  about  orthogonal 
transformations.  Every  step  except  two  and  four  either  does  not  affect  F  and  G  or  just 
multiplies  them  on  the  left  or  right  by  orthogonal  transformations,  thus  for  those  steps 
we  need  only  invoke  orthogonal  invariance. 

For  Step  2,  we  have  that 


r*(2)  __  j?(l) 


-  e  eT  F 

rr 


(i) 


\5n  / 


=  ui 


\  9n  / 


Thus  ||(F^  G<2))||f  <  ||(F^  follows  from  the  fact  that  F^  and  G ^  are 

no  larger  than  F^  and  componentwise.  For  the  Euclidean  norm,  we  note  that  for 
any  vector 


v  — 


k 


p-k 


52 


thus 


\\(F{1)  G(1)H |2  =  ||(-P(2)  G(2))^||2  +  ({/(1)}:rf(1))2 

+  ({^V  -  {g®}2)(  v(2b2 

>  ||(-P(2)  G(2))v||2 


for  all  v  G  7 Zn.  A  simple  argument  on  the  definition  of  the  Euclidean  norm  yields  the 
inequality  ||(F^  G^)||  >  ||(^2^  G^)||. 

Step  4  clearly  produces  |<7^|  <  The  same  argument  as  above  yields 


||(F(2)  G(3))||  <  IK^2)  G(2))|| 


and  likewise  for  the  Frobenius  norm.  Thus  we  have  (4.4).  □ 

In  the  absence  of  rounding  error,  we  can  show  that  the  “additive  noise”  in  Step  4 
satisfies  an  important  consistency  property.  We  assume  that  C  is  orthogonally  equivalent 
to  A  +  SA  where  6  A  consists  of  errors  from  previous  orthogonal  transformations. 

Proposition  4.2.  Assume  that  Algorithm  J^.l  is  performed  in  exact  arithmetic ,  that  U 
and  V  in  (4-1)  are  exactly  orthogonal.  Let  U  and  V  be  as  in  (1.13)  and  let  z  =  V  r  be 
as  in  (4-7)-  Assume  that  Ue^  =  e  ,  and  that  z  =  V  r  is  computed  exactly.  Also  assume 
that  C  satisfies  (4-1)  with  backward  error  SA,  that  is 


A  +  SA  = 


T  £  T  \ 
r  +  or 

=  u 

\ 

c 

^  A  +  SA  j 

V  0  ) 

(4.8) 
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Then  z  and  C  satisfy 


(4.9) 


Thus  ||M  ||  =  ||«A||  <  ||M||  and  the  same  result  holds  for  the  Frobenius  norm . 
Proof.  We  have  that 

T 

x  -f  6x  ( p  —  8p)e  0 

he ?  o 

jp(2)  0 

0  0  0 

We  also  have  that 


c  ^ 

0  ldiae(4'vi-W  = 
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Let  Vj  =  V diag(/fc,  )  then  from  the  definition  of  V  in  Algorithm  4.1  it  follows 

that 


(4.10) 


However,  we  have 


Comparing  equations  (4.10)  and  (4.11),  using  the  fact  that 

x 

p6l 

yo 

and  the  assumption  that  Ue  =  e^,  we  have  (4.9).  □ 

Therefore,  we  have  shown  that  the  additive  noise  in  Step  4  and  the  act  of  ignoring 
2/q  actually  make  the  matrix  C  closer  to  being  orthogonally  equivalent  to  A  than  C  is  to 
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A.  In  Section  4.3,  we  show  that  the  results  of  Propositions  4.1  and  4.2  make  Algorithm 
4.1  a  robust  algorithm  for  tracking  the  ULVD. 

4.2.2  Relation  to  Park  and  Elden’s  URV  Procedure 

A  recent  report  by  Park  and  Elden  [87]  gives  a  method  for  downdating  the  URVD. 
For  that  algorithm,  we  are  downdating  an  upper  triangular  matrix  M  of  the  form  (1.4). 
The  procedure  is  as  follows. 

1.  Find  an  orthogonal  matrix  U  and  an  upper  triangular  matrix  R  such  that 


(*T\ 


\ 


=  U, 


R 


\ 


Determine  S  and  y  such  that 


(vT\ 


(4.13) 


(4.14) 


2.  Find  an  orthogonal  matrix  U  and  an  upper  triangular  matrix  T  such  that 


_T  \ 

\ 

y 

0 

=  02 

V  f  ) 

The  downdated  matrix  M  is 


I  R 


(4.15) 
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as  before. 

Park  and  Elden  [87]  recommend  the  use  of  hyperbolic  rotations  in  (4.14).  That 

rp  rp 

can  be  avoided  by  a  simple  and  well-known  trick.  Let  (a^  )  be  the  first  column  of 

U  as  determined  in  (4.13).  Then  we  note  that 


/  \ 


y  =  (y  s  ) 


O' 


which  implies  that 


y  =  a1  l(y-STal). 


Once  y  is  determined,  we  obtain  S  from 


I  T\ 
\  S  } 


=  u. 


t  f\ 

\  5  / 


T  T 

If  T  T  -  yy  is  indefinite,  Park  and  Elden  substitute  for  (4.15)  the  operation 


N 

(  T  \ 

0 

pe 

V  = 

n 

f  j 

_  71  _ 

where  V  is  an  orthogonal  matrix  such  that  V  y  =  and  so  that  T  is  triangular.  They 
then  compute 


T  \ 

T  \ 

y 

t-tT 

y 

=  U1 

\ s  > 

\ s  > 
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where  y  =  V  y.  The  downdated  matrix  T  =  T  except  for  t  ,  q  =  n  -  k.  That  entry 


satisfies 


f  1 0 _  if  i‘„i  <  w 

^  |  ]j^qq  ~  ^  ot^6rw'se  • 


j1  _ t  t 

The  first  condition  is  identical  to  the  case  where  T  T  —  yy  is  indefinite. 


Remark  4.3.  Algorithm  4.1  is  related  to  the  above  algorithm  although  the  essential 
steps  of  them  were  developed  independently.  Steps  1  and  2  of  Algorithm  4.1  reduce  the 
ULVD  downdating  problem  to  that  of  downdating  the  (k+  1)  x  (k  +  1)  upper  triangular 
matrix 


k  l 


L ^  h  |  * 


•  *s\. 


(4.16) 


The  Park-Elden  algorithm  applied  to  this  matrix  would  perform  Steps  3  and  4  of  Algo¬ 
rithm  4.1  except  that  it  would  perform  Step  4  according  to 


yT  I  P 


(4.17) 


instead  of  according  to  (4.6).  The  consistency  property  in  Proposition  4.2  does  not  hold 
if  we  use  (4.17).  That  can  be  seen  when  (4.17)  is  placed  into  the  proof.  Equation  (4.10) 
does  not  change,  but  we  do  have 


p  - 

h  —  [6p\a 


(4.18) 


_  _  ^2)  —1 

where  Sp  =  p  -  '  =  a  tip  using  terminology  from  Algorithm  4.1.  We  note  that  the 

additive  noise  is  multiplied  by  a  factor  of  a  *  >  1,  but  that  is  not  crucial.  Equation 
(4.18)  implies  that  equation  (4.11)  in  the  proof  of  Proposition  4.2  is  replaced  by 

x  +  Sx  (p  —  [6p\a)e^  0 

(h  —  [Sp\a)e^  0 
^(2)  g(3)  0 

0  0  0 

Thus  the  consistency  property  of  Proposition  4.2  does  not  hold. 

The  Park-Elden  URVD  algorithm  requires  5fc2  +  5(n-fc)2  +  8fc(n  —  k)  +  0(n)  flops. 
This  algorithm  is  fewer  flops  than  is  required  to  downdate  the  ULVD  by  Algorithm  4.1, 
but  there  is  an  important  advantage  to  maintaining  the  ULVD  instead  of  the  URVD. 
Both  the  URVD  and  the  ULVD  will  tend  to  produce  a  V  matrix  such  that 


V’  vie 


v2 «  K 


where  range(Vj)  and  range^)  are  approximations  to  the  subspaces  associated  with  the 
first  k  and  last  n  -  k  singular  values  of  A  respectively.  However,  as  found  by  Fierro  and 
Hansen  [43],  if  F  in  (4.2)  and  S  in  (1.4)  satisfy  ||F||  &  ||5||  and  if 


inf  ll^ll  >  e  >  .max  ||Gy||, 

11*11=1  lly||=1 


then  the  ULVD  yields  a  more  accurate  approximation  of  the  desired  subspaces  of  A . 
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4.3  Error  Analysis 

4.3.1  Error  Bounds  on  Algorithm  4.1 

~  (2) 

We  begin  by  showing  that  if  p  <  gKu'  in  Step  4  and  yQ  =  0,  then  Algorithm  4.1 
produces  a  matrix  C  such  that 

VTCTCV  ~(z  +  6z)(z  +  Szf  =  (C  +  SC)T(C  +  SC)  (4.19) 

for  some  orthogonal  matrix  V.  This  is  the  so-called  mixed  stability  criterion  defined  in 

(2) 

Definition  2.7.  If  in  Algorithm  4.1,  yQ  ^  0  or  p  >  then  we  are,  in  fact,  producing 
C  such  that 

VTCTCV  ~(z  +  Sz  +  SzQ)(b  +  6z  +  SzQf  =  (C  +  SC)T(C  +  SC) 
for  some  orthogonal  matrix  V.  In  the  context  of  Algorithm  4.1  , 


We  note  that  V  =  I  ,  if  Step  6  is  not  done.  The  effect  of  Sz  is  discussed  in  Appendix 

o  p  1C  U 

B  of  [13].  It  is  essentially  the  same  as  the  effect  of  6b  as  bounded  in  Proposition  4.3. 
The  analysis  in  this  section  will  ignore  that  effect,  it  will  assume  that  we  are  analyzing 
the  problem  (4.19). 
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where  0^(p)  is  a  modestly  growing  function  of  p. 

In  the  next  section,  we  discuss  the  effect  of  these  rounding  errors  on  the  singular 
vectors  of  the  matrix  C. 

4.3.2  Effect  of  Rounding  Errors  on  Singular  Vectors 

A  common  reason  for  computing  the  ULVD  is  to  separate  subspaces  associated 

with  large  and  small  singular  values.  For  the  ULVD,  after  downdating  ,  there  will  be  l 

large  singular  values  where  /  =  k  or  /  =  k  —  1.  If  W  —  (W^  W^),  W ^  £  1ZnX\  W ^  £ 

TZ nx(n-l )  .g  tjie  matrjx  Df  j-jght  singular  vectors  of  C,  then  the  computed  range(W1) 

and  range(W  )  should  be  as  close  as  possible  to  the  expected  ranges.  In  this  section, 

we  show  how  reliable  we  can  expect  these  subspaces  to  be.  Our  bounds  are  somewhat 

better  if  /  =  k,  that  is,  if  the  downdate  does  not  alter  the  rank. 

We  need  to  measure  the  effect  of  both  6  and  6  on  the  invariant  subspaces 

z  C 

associated  with  the  l  largest  singular  vectors  and  the  n  —  l  smallest  ones.  We  define  C 
and  C  as  matrices  such  that 


(4.23) 

(4.24) 

Defining  C  as  in  (4.19)  implies  that 


CTC  =  VTCTCV  -  zz 


CTC  =  VTCTCV  -  (z  +  Sz )(z  +  6z)T. 


C  =  C  +  6C. 


(4.25) 


It  is  also  necessary  to  define  u  by 


u  =  max{l,||[I(2)]-Tx||,||i-:rx||}. 


(4.26) 


Note  that  the  definition  of  u  involves  only  the  L  block  and  is  related  to  the  condition 
number  given  by  Pan  [86]. 

First,  we  need  the  following  three  technical  lemmata. 

Lemma  4.1.  The  vector  h  and  scalar  p  from  Algorithm  4-1  satisfy 

<  Ill'll-  (4.27) 


M 


Proof.  We  note  that 

/  \ 

P 

\  1 

which  completes  the  proof.  □ 

Lemma  4.2.  Let  z  be  as  in  (4-7)  resulting  from  Algorithm  4-1 ■  Then 

IM<l|G||(l+u). 

Proof.  From  Algorithm  4.1  and  (4.20)  we  have  that 


\ 

\ 

X 

X 

Z  — 

V 

=  diag(4'  vl  •  W  di^(«f  ■  • 4-i-i  >  diag(4'  •  4 -,) 

y 

<  4+i  °  >pT 


'o' 


V. 


<  INI, 
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which  completes  the  proof.  □ 

It  is  necessary  to  bound  the  effect  of  A  and  z  on  the  singular  vectors  of  C. 

LEMMA  4.3.  Let  w.,  i  =  1,2 be  the  right  singular  vectors  of  C .  Let  A  be  defined 
by  (4-21),  let  z  be  as  in  Lemma  4-2,  and  let  u>  be  given  by  (4-26).  Then 

l|A».||  <  l|X_I||  \/»2  +  l|G||2<l|i‘1||  (ff.  +  IIGII),  i=  1,2,.. 

\zTw.\  <u  (ff.  +  2  ||G||),  i  —  1,2, . .  .,p. 


Proof.  We  note  that  for  i  =  1, 2, . . . ,  p, 


It.is  obvious  from  the  form  of  C ,  that  the  last  n  —  p  components  of  all  its  nonzero  singular 
vectors  will  be  zero.  It  should  also  be  noted  that  w.  =  e.,  i  =  p+ 1, . . . ,  n  form  a  singular 
subspace  for  the  last  n  —  p  zero  singular  values.  For  i  =  1, 2, . . . ,  p,  we  have  that 


||A«,.||2  =  ||K<1>||2  +  ||G||2||i-1||2||»<2)||2. 


Taking  square  roots  establishes  the  first  bound  on  ||Aw.||.  Continuing  we  have 

||At».||2  =  ||i-1Z«.<1)||2  +  IIG|l2lli'1ll2 
<  ||i'1||2(||Zti>|I)||2  +  ||G||2) 
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<  P“1II2(^  +  II0||2) 

Taking  square  roots  again  establishes  the  second  bound  on  ||Ato.||. 
Now 


,-T  ,  .  I -T  (1),  ,  ,_T  (2), 

\z  w.\  <  \x  w)  |  +  |t/  w.  '| 

<  |xTX-1Z^1)|  +  \\y\\ 

<  iirtii  <r. + m 

<  u.  (<r.  +  ||G||)  +  lies'll 

<  u>  (a.  +  2  lies'll), 

which  completes  the  proof.  □ 

Now  we  give  a  perturbation  bound  on  the  effects  of  the  backward  errors  ||<5z||  and 

Ml- 

Lemma  4.4.  Assume  the  terminology  of  Lemma  f.3.  Let  w w.,  w i  =  1,2, ...,n  be 
the  right  singular  vectors  of  C,  C,  and  C  in  (4-19),  (4-%3)>  and  (4-%4)  respectively.  Let 
cf  j  >  •  •  •  >  <Tn  be  the  singular  values  of  C ,  and  let  >  •  •  •  >  a ^  be  the  singular  values 
of  C.  Ifo.yo.  and  o.  >  a .  we  have 


\w.  W.\<u  \\Sz\\ 


+  <?(M2) 


\wTw.\  < 
1  i  3 


+  0(\\sc\\2). 


(4.28) 


a .  —  a . 

i  j 


(4.29) 
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Proof.  Equation  (4.29)  is  just  a  standard  error  bound  based  upon  the  perturbation 


z=,T 


of  the  eigenvectors  of  C  C.  From  the  Kato  [69]  expansion  for  eigenvectors  we  obtain 


\vk  w  ,\  = 

*  r 


< 


wT  (z(6z)T  +  (Sz)z^  +  (6z)(Sz)^)w . 
_ l _ J 


~2  ~2 
a.  —  a . 
*  J 


+  °(INI2) 


INI  (lf'r’5il  +  lfT'5'l) 


~2  ~2 

a .  -  a . 
i  3 


+  0(\\6zf). 


From  Lemma  4.3  we  have 


\~T , 


a  +  a  +  4  ||G|| 


|u>tf  u?  J  <  11^1  1  _2J  _2  -  +  0(\\6z\f). 


a .  —  a . 
*  3 


An  algebraic  simplification  leads  to  (4.28).  □ 

Simply  using  the  definition  of  the  ||  •  || ^_,  norm  leads  to  the  following  proposition. 
It  tells  us  how  good  our  subspaces  will  be  if  we  only  have  a  bound  of  the  form  (4.22). 

PROPOSITION  4.3.  Assume  the  terminology  of  Lemma  4-3.  Let  W  =  ( W  W  ),  W  E 

12  1 

Wy  E  7Znx(n  0  be  the  matrix  of  right  singular  vectors  of  C  and  let  W  = 
(W^  Wy)  be  the  corresponding  matrix  of  right  singular  vectors  of  C,  If  c >  <7j 
then 


ll«f  W2\\F  <  n  -  I) 


U  1IM1  +  Ill'll 

^  ~  ai+ 1 


4W||G||  ||*S|| 


-1 

\ 

2\ 

Sz 

+  o 

[  SC  ) 

/ 

Proposition  4.3  applies  even  if  the  downdate  changes  the  rank.  The  condition 
number  u  in  (4.26)  depends  solely  upon  the  L-block,  the  matrix  C  does  not  have  to  be 
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well-conditioned  or  even  full  rank  for  the  downdating  problem  to  be  well-conditioned  as 
is  required  in  previous  analyses  [86,  87].  If  there  is  no  rank  change,  that  is  if  l  =  k,  we 
can  get  an  even  better  bound  as  shown  below. 

We  define  e  and  E  according  to 

e  =  A ~lSz/6z,  E  =  6C  (4.30) 

where  6  and  are  the  bounds  from  Theorem  4.1.  Note  that  ||e||  =  ||£||  =  1.  First,  a 
z  O  t 

technical  lemma  characterizes  the  vector  z  that  is  downdated. 

LEMMA  4.5.  Assume  the  results  and  terminology  from  Theorem  4-1-  Let  uL,  i  =  1, 2, . . .,  n 

be  the  right  singular  vectors  of  C  and  let  t».  ,  i  =  1, 2 ,...,«  be  the  right  singular  vectors 

of  C.  Let  a  >  cf2  >  . . .  >  be  the  singular  values  of  C.  Then  for  all  i  and  j  such  that 

a .  ±  c .  we  have 
i  3 

I w  -I  <  _2  Z„2(IIA%-11  +  llA55ill  +  °(^)-  (4-31) 

i  3 


Proof.  From  Kato[69],  we  have  that 


~T  ~ 
w .  w  .  — 
*  3 


~T  _  T  -T  ~ 
w  .  A (ze  +  ez  )Aw . 

3  i 


+  o(s2z) 


where  e  is  defined  in  (4.30).  Taking  norms  we  obtain  (4.31).  □ 

Lemma  4.6.  Assume  the  results  and  terminology  from  Theorem  4-1.  Let  w.,  i  =  1, 2, . . .,  n 
be  the  right  singular  vectors  of  C  and  let  w.  ,  i  =  1,2, . .  .,n  be  the  right  singular  vectors 


\ 


Reorganizing  leads  to  (4.32).  □ 
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Now  we  can  bound  the  quantities  in  our  lemma  on  the  singular  vectors.  These 
bounds  actually  apply  to  the  singular  vectors  of  C,  but  we  will  ignore  second  order  effects 
and  use  these  bounds  for  the  singular  vectors  of  C. 

Now  we  need  a  bound  on  the  errors  in  the  singular  vectors  of  C  from  6  .  To 
within  rounding  error,  the  same  bounds  as  in  Lemma  4.3  will  hold  for  uL,  i  —  1, 2, . . . ,  n. 
That  is,  \\zTw.\\  <  u  (a.  +  2  ||G||)  and  ||A5L||  <  \\L  *11  (a.  +  ||G||). 

Lemma  4.7.  Assume  the  terminology  of  Lemma  J^.3  and  the  results  and  terminology  of 

Theorem  A.L  Assume  that  a.  >  a ..  Then 
*  1  3 

S  lo  II L  ^11 

I  <0,1  <2  ‘  ■"  5  +  1.5  ||G||  +  2  ||G||2/(5  +  J )  +  0(i\ 

J  1  O .  —  O  .  J  J  1  6 


Proof.  Combining  Lemmata  4.3  and  4.5  yields 


r  6  w  II1  *11  r 

wjw.\  <  *  2  2  [(gt-  +  2  Hgll)(gJ-  +  llGll) 

°i  ° j 

+  (a.  +  2  IIGIIXSj  +  IIGII)]  +  0(1]) 

5  UJ  ||  X I  j| 

=  2  2  2  (Vi  +  1-5  ||G||(ff.  +  S  )  +  2  ||G||2)  +  0<) 

0,-0.  J  J 
1  3 

6  u>  ||X  1  j]  r  o  1  2 

<  2  -  <7  +  1.5  ||G||  +  2  ||G||2/(cr  +  cr  )  +0(62), 

o .  —  o .  3  1  J  z 

i  3 


which  completes  the  proof.  □ 

The  proof  of  the  following  theorem  is  obvious  from  Lemma  4.7. 
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THEOREM  4.2.  Assume  the  hypothesis  and  terminology  of  Lemma  4.7.  Let  W  be  the 
matrix  of  singular  vectors  of  C  and  let  W  be  the  matrix  of  singular  vectors  of  C.  If 


W  =  (WX  W2),  W  =  (W1  W2), 


where  W  %W.  G  TZnxk,  and  W  ,W  G  Knx(n  k )  then 


_ i -  8  LO  ||i  II 

\\<W\\  <  2fkin-k)^—= - X 

0k  k  + 1 

h+ 1  +  1'5  l|0|l  +  2||G||2/(5t  +  5t+1)]+0(«22) 


Analogously,  we  can  bound  the  effect  of  <5^. 

Theorem  4.3.  Assume  the  hypothesis  and  terminology  of  Lemma  J^.6.  Let  W  be  the 
matrix  of  singular  vectors  of  C  =  C  +  SC  and  let  W  be  the  matrix  of  singular  vectors  of 
C.  If 

w  =  (Wl  w2),  W  =  (Wl  W2), 

where  W  W  G  TZnxk ,  and  W  W  G  then 

11  z  z 


II  wTw. 


2  "F 


<  2  fk{n  —  ~k) 


II  (* 


fc+1 


ak  ak+l 


+  I|G||)  o 

- +  o(s2c). 


The  effect  of  6 _  is.  in  fact,  somewhat  less  critical  than  that  of  S  as  has  been 
C  z 

stated  in  other  analyses  of  this  problem  [86,  87,  99].  We  note  that  the  error  bounds  in 
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Theorems  4.2  and  4.3  are  relative  gap  bounds  on  the  error  in  the  subspaces  similar  to 
those  in  [11,  36]. 

If  ||Z-1||  ||G||  <  1,  these  bounds  are  a  significant  improvement  over  those  in 
Proposition  4.3.  This  is  one  of  the  reasons  for  maintaining  the  property  (4.4). 

4.4  Numerical  Examples 

In  this  section,  we  present  a  few  examples  from  numerical  experiments.  These 
tests  were  performed  using  Matlab  on  a  SPARCstation  5  in  IEEE  Standard  double 
precision  with  machine  precision  «  10  .  The  algorithm  employed  the  sliding  window 

technique  described  in  Section  2.7.3. 

At  each  step  of  the  sliding  window  method  with  the  window  size  mQ,  an  X  n 

data  matrix  is  constructed  from  an  m  X  n  observation  matrix  A  by  adding  a  new  row  to 

the  data  matrix  in  the  previous  window  and  deleting  the  oldest  row  from  it.  In  step  j, 

the  row  +  j  of  the  observation  matrix  is  added  and  the  row  j  is  deleted,  giving  the 

data  matrix  A .. 

J 

Computing  the  ULVD  of  initial  data  matrix  is  described  in  Section  2.5.2.  Then 

Algorithm  4.1  takes  the  lower  triangular  matrix  (middle  part  of  the  decomposition), 

the  orthogonal  matrix  (right  part)  as  initial  input  and  the  modifying  vector  r,  and 

T 

successively  modifies  these  matrices  at  every  window  step.  The  vector  z  =  V  r  is 
computed  at  the  beginning  of  each  window  step. 

We  tested  our  algorithms  in  the  context  of  the  total  least  squares  (TLS)  problems. 
See  Section  2.7  for  details.  We  used  the  TLS  solutions  from  the  Jacobi  SVD  as  reference 
in  checking  the  accuracy  of  the  solution  and  rank  estimates  of  our  algorithms. 
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Fig.  4. 3-4. 5  show  the  rank  estimates  by  our  algorithm,  which  are  identical  with 
those  of  the  Jacobi  SVD  algorithm.  The  horizontal  axis  represents  the  window  steps  and 
the  vertical  axis  the  numerical  rank  of  the  window  matrix. 

The  distance  between  the  subspaces  is  given  in  the  next  plot  using  the  Definition 
2.3.  Let 

A.  =  Y.X.WT ,  W .  =  (W.  W.n) 

3  3  3  3  3  v  3 2' 

be  the  SVD  of  A.  computed  by  the  one-sided  Jacobi  method  at  step  j.  Let 

A.  =  U.C.VT ,  V.  =  (V..  V.A 
3  3  3  3  j 

be  the  ULVD  of  A.  computed  by  Algorithm  4.1.  Note  that  here  we  are  discarding  Uy 
Finally,  let 

C.  =  ?.Z.WT,  W.  =  (W.,  W.J 
3  3  3  3  3  j2 

by  the  SVD  of  C .  computed  by  the  one-sided  Jacobi  method.  Define  W .  by 

W.  =  (W.,  W.\  =  V.W.. 

3  Jl  3%  33 

Define  the  angles  between  the  subspaces 

=  sin(»2)  =  ||^F.2||,  Sm(»3)  =  (4.34) 

The  angles  6.,  i  =  1,2,3  represent,  respectively,  error  between  the  true  noise  sub¬ 
space  from  the  Jacobi  SVD  and  the  approximate  one  from  tracking  the  ULVD,  the 
approximation  error  from  the  ULVD,  and  the  subspace  errors  from  ULVD  subspace 
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tracking.  The  value  sin(^)  is  one  that  is  bounded  by  our  error  analysis.  We  plot¬ 
ted  log10(sin(6h)),i  =  1,2  in  solid  and  dotted  lines,  respectively,  in  the  vertical  axis  of 
the  second  graph.  It  turned  out  that  log10(sin((?3))  was  almost  indistinguishable  from 
loglo(sin(01)),  so  we  did  not  plot  it.  sin(02)  is  the  approximation  error  discussed  by 
Fierro  and  Hansen  [43]. 

Finally,  the  TLS  errors 


||I(SVD)_I(UIVD)|| 

2  life’ll 

are  given  in  logarithm  in  the  last  plot.  Here,  x and  are  the  TLS  solutions 

using  the  SVD  and  the  ULVD,  respectively. 

For  our  condition  estimators,  we  use  the  LINPACK  condition  estimator  to  approx¬ 
imate  the  left  singular  vector  that  corresponds  to  the  smallest  singular  value,  followed  by 
inverse  iteration  using  this  approximate  singular  vector  as  the  initial  vector.  The  tests 
show  that  up  to  three  steps  of  inverse  iterations  suffice  the  accuracy  of  the  approximate 
smallest  singular  value  required  by  the  algorithm. 

Example  4.1.  A,  a  110-by-6  random  matrix,  6,  a  110-by-l  random  vector.  Entries  of  A 

and  b  were  chosen  from  a  uniform  distribution  on  the  interval  (0, 1).  85  randomly  chosen 

rows  of  (A;  6)  were  multiplied  by  7  =  10  in  order  to  vary  the  rank  of  the  matrix,  and 
—2 

tol  =  10  .  The  window  size  p  used  was  12. 

_ g  _ 0 

EXAMPLE  4.2.  Same  as  Example  4.1  except  that  7  =  10  and  tol  =  10 

4 

EXAMPLE  4.3.  Same  as  Example  4.1  except  that  the  matrix  had  an  outlier  of  size  10 
at  (18, 1)  position. 
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The  first  plot  shows  that  our  algorithm  estimated  the  numerical  ranks  correctly 
throughout  the  sliding  window  steps  in  spite  of  frequent  rank  changes.  Although  the 
errors  in  tiny  singular  values  were  relatively  large,  and  the  small  singular  values  were 
almost  always  overestimated,  they  were  close  enough  to  correctly  estimate  the  rank.  Thus 
the  rank  estimates  from  our  algorithm  and  those  by  the  Jacobi  SVD  were  identical.  As 
expected,  the  errors  in  the  TLS  solution  t.  are  almost  exactly  the  same  as  the  size  of 
sin(0j)  in  (4.34). 

The  second  plot  in  each  figure  shows  that  the  noise  subspace  error  is  very  small 
giving  accurate  TLS  solutions.  The  quantities  in  (4.34)  are  shown  to  be  essentially 
identical  indicating  that  the  subspace  errors  from  our  algorithm  are  from  the  rounding 
errors,  not  approximation  errors.  We  note  that  the  Example  4.1  has  greater  error  in  the 
noise  subspace  than  Example  4.2. 

This  is  probably  because  Example  4.1  has  only  a  small  relative  gap  in  the  spectrum 

—5  _c 

around  7  =  10  but  a  large  relative  gap  around  7  =  10  .  However,  for  all  of  our 

examples,  the  approximation  to  the  subspaces  by  the  ULVD  is  very  good.  Since  the  error 

bounds  on  the  distance  between  the  noise  subspaces  depend  on  the  ( k  +  l)-st  singular 

value,  the  approximated  singular  subspace  gets  better  as  7  decreases  as  shown  in  Fig. 

4. 3-4. 4. 

The  TLS  solution  errors  behave  very  much  the  same  as  the  noise  subspace  errors. 
From  (2.18)-(2.19)  it  is  not  difficult  to  see  that  the  TLS  errors  differ  from  the  noise 
subspace  errors  only  by  a  constant  factor. 

Moreover,  the  algorithm  performs  well  even  when  G  becomes  singular  (indicated 
by  in  the  first  plot).  We  tested  several  other  examples,  and  these  results  were  typical. 

Since  our  downdating  procedures  use  the  LINPACK  downdating  algorithm,  it  is 
not  difficult  to  generate  the  cases  where  the  algorithm  breaks  down  when  ||a||  >  1,  for 
instance,  when  deleting  a  row  that  contains  outliers.  In  this  case,  we  first  refine  the 
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decomposition  [100,  104],  that  is,  compute  an  orthogonal  matrix  U  such  that 


~  ~  \ 

L  0  1 

L  H 

u 

= 

^  G  ) 

K°  G  ) 

(4.35) 


~T 

If  it  is  still  true  that  ||a||  >  1  even  after  the  refinement,  where  L  a  =  x,  the 
corrected  semi-normal  equation  (CSNE)  approach  [18]  (indicated  by  ’+’  in  the  first 
plot)  is  used  for  computing  a  with  higher  accuracy.  It  is  essentially  the  same  as  that 
used  by  Park  and  Elden  [87]  and  is  given  by 


Ly  =  a,  t  =  e1-  X  .V^y 


I?  8a  =  vJxTt,  a  =  a  +  8a 
1  J 

L8y  =  8a,  t  =  t  —  X  ,V^8y,  a  =  ||f|| 


where  xT  =  aJj  ,  the  j-th  window  matrix  augmented  with  the  row  being  deleted. 
Finally,  restore  the  lower  triangular  form,  that  is,  compute  an  orthogonal  V  such  that 


- 

\ 

~ 

~  \ 

L 

0 

L 

H 

d) 

1° 

G  ) 

(4.36) 


The  CSNE  approach  was  used  in  all  three  examples  and  most  extensively  in 

Example  4.3  when  downdating  a  row  with  an  outlier.  However,  the  performance  of 

our  algorithm  was  less  satisfactory  for  larger  outliers.  This  is  consistent  with  our  error 

analysis,  since  for  a  very  large  outlier  we  would  have  u  in  (4.26)  very  large. 

The  refinement  steps  in  (4.35)-(4.36)  require  k  (n  -  k)  Givens  rotations,  so  that 

1 

they  become  impractical  when  k  =  0(n2)  or  larger.  As  an  alternative,  we  solve  for  a 
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Table  4.1.  Tracking  || jP||  and  ||G||  for  the  ULVD  Procedure 

F  r 
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Moreover,  to  prevent  a  from  becoming  too  large,  we  track  ||P||^,  so  that  it  remains 

~1 

under  certain  threshold,  say,  \\L  ||  ||.F||  <  0.01.  This  is  similar  to  recommendations 

made  by  Fierro  and  Bunch  [41,  42].  Only  steps  in  Algorithm  4.1  that  require  to  update 

||JF||^,  are  Steps  2  and  5.  Table  4.1  shows  how  to  update  these  quantities.  Here, 

(i) 

denotes  the  first  column  of  GK  ' ,i  =  1,3,4. 


*** 


Rank  Estimates 
**  ***** ***** *  **  '  **7 

+  + 


Noise  Space  Errors 


t - 1 - 1 - r 


U  •  *N  . 


I _ I _ I _ I - 1 - 

0  30  40  50  60  7 

TLS  Solution  Errors 


i 


30 


40 


50 


60 


Fig.  4.3.  Example  4.1 


Fig.  4.5.  Example  4.3 
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Chapter  5 

Rank  Detection  for  Modifying  the  ULV  Decomposition 

5.1  Introduction 

In  Algorithm  4.1  the  deflation  steps  are  required  to  compute  the  numerical  rank 
after  a  downdate.  The  deflation  step  usually  involves  some  condition  estimator,  which 
can  be  a  nuisance  in  some  situations  such  as  in  parallel  implementation.  A  survey  of 
popular  condition  estimators  is  given  in  [61],  and  they  all  require  0(k  )  flops,  so  that 
the  entire  process  requires  0(k  )  flops,  where  k  is  the  dimension  of  L  in  (4.2).  With 
some  modification  to  Algorithm  4.1  we  can  often  eliminate  Step  7,  the  deflation  step. 

Furthermore,  the  modified  algorithm  offers  an  efficient  way  of  tracking  exact  quan¬ 
tities  of  || L  *||^,  IlF’H^,  and  ||G||^,,  which  give  a  significant  information  on  the  condition 
of  the  downdating  problem.  The  experiments  show  that  the  computed  subspaces  are  as 
good  as  can  be  expected,  and  no  worse  than  those  demonstrated  in  the  previous  chapter. 
The  updating  algorithm  for  the  ULVD  can  be  also  implemented  similarly  with  a  slight 
modification  to  the  downdating  algorithm. 

We  propose  our  new  ULVD  updating/downdating  algorithm  in  Section  5.2.  Sec¬ 
tion  5.3  contains  the  rank  detection  method  related  to  the  new  algorithm.  In  Section 
5.4,  we  give  numerical  tests  of  our  algorithm  on  the  RTLS  problems. 
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5.2  New  Algorithm  for  ULVD  Downdating 

Our  new  algorithm  reduces  the  downdating  problem  to  a  2  x  2  downdating  prob- 

T  T 

lem.  As  in  Algorithm  4.1,  we  assume  that  L  L  -  xx  is  positive  definite,  and  U  is  not 
accumulated. 

Algorithm  5.1  (New  Procedure  for  ULVD  downdating).  Assume  the  terminol¬ 
ogy  of  Algorithm  4.1,  and  denote  also  that 

Step  1.  Construct  orthogonal  matrices  €  TZ^X^  and  U^,  V,  £  ^)x(p  ^) 

such  that 


Also,  compute 


(1). 
t  - 

-DxLVv 

JL 

II 

(5.1) 

(1)  = 

-u^av2, 

V*y  =  pev 

P  =  Ill/ll- 

(5.2) 

F(1)  =  U^FVy 


Update  V  *-  V  d iag(Vj,  J^)  diag^,  Vyl  ). 


Step  2.  Find  an  orthogonal  matrix  U  G  7jd"bl)xd+l)  such  that 


L^>  h 


0  9] 


.T  [  i(1)  0 

73  rf(DyT  JD 
\  *  9n 


Define  P(2)  =  (I  -  . 
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Step  3.  Write 


r(2)_ 


ifc-1  1 

l  £®  0 

,  (2)tT  ,(2) 

>  hi 


fc— 1 


jkxjfc 

and  find  an  orthogonal  matrix  such  that 


£®  = 


/  £®  «,<3>  \ 


0 


Z(3)  ; 

kk  } 


4 


(5.4) 


m  -T 

Compute  also  hK  '  =  1/  h . 


Step  4.  Use  Algorithm  3.9  for  2  x  2  downdating: 


where  we  solve 


/  (i3)  W  f  N  /  - ' 


0  5 


(2)  / 
11  / 


'1 

V  ^2  / 


* 


(5.5) 


and  set 


P \~P^  P2~  min^2’  ^ 


(5.6) 


kxk 

Step  5.  Find  an  orthogonal  matrix  V  €  72-  such  that 


z(5)  = 


r(5) 

^ll 
r  (5),T  ,(5) 

'}  i, 


=  l<4V„ 


jfcjfc  / 


(5.7) 
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where  L ^  and  L ^  differ  only  on  the  (k,k)  entry.  Also,  compute 


p(3)  _  p(fyy 

3‘ 


Update  V  «-  Vdmg(V3,In_k). 


Step  6.  Find  an  orthogonal  matrix  V  €  7£(k+l)x(k+l)  such  that 

(T(6)  0)  =  (I(5)  hi2)e^)V4 
(F{4)  G(4))  =  (F^  C?<3>) diae(^4, J  *_!>- 


Update  V  «-  Udiag(U4,/n_]b_1).  If  ±  0,  then 


10  0 


0  0 


c  =  F  G  0  =  G^  0 


0  0  0 


0  0  0 


and  skip  Step  7. 


Step  7.  If  =  0  then  =  0  also  since  it  was  formed  from  using  V^.  Thus 


k  p-k 


k  0 


p-k  I  j?(4)  (j(4) 


k  p—k 


L  o 


=  I  o  o 


F  /  P-k-i 
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~(4) 

where  GK  '  is  a  lower  Hessenberg  matrix.  We  then  find  an  orthogonal  matrix 
V  £  7Z^P  k)*(p  k)  an(j  an  orthogonal  permutation  matrix  U  £  Tl^P-k)x(p-k) 

o  0 


such  that 


°Ig{\= 


<6  o' 

l"'/ 


p—k—1 

1 


Update  V  «-  V  diag (/, ,  V_,  I  ).  Thus, 

K  O  Tl — p 


k  p-k-1  n-p+1 


>i 


c  = 


0  0 

F  G  0 

0  0  0 


k 

p—k— 1 
n— P+1 


Step  8.  Update  the  numerical  rank  of  L,  and  make  appropriate  adjustments  to  F 
and  G  (a  procedure  for  the  rank  detection  will  be  described  in  the  next  section). 

Algorithm  5.1  requires  12 fc2  +  6(p  —  k)^  +  24 k(p  -  k)  +  0(p)  flops  for  Steps  1-6  and 
4k  +  0{k )  flops  for  Step  7.  When  V  is  modified,  additional  6 np  +  12 nk  flops  will  be 
required.  Fig.  5. 1-5.2  depict  the  reduction  steps  for  Algorithm  5.1.  As  in  Algorithm 
4.1,  a  set  of  —>  denotes  the  downdating  Step  4,  and  Step  7  is  illustrated  in  Fig.  4.2. 

Remark  5.1.  Note  that  in  (5.5)-(5.6)  we  incorporate  “additive  noise”  to  compute  a  as 
similarly  done  in  Step  4  of  Algorithm  4.1.  As  shown  in  Chapter  4,  the  algorithm  can  be 
made  consistent  with  this  additive  noise  in  the  absence  of  rounding  error.  Proposition 
4.2  shows  this  important  consistency  property. 

Remark  5.2.  Algorithm  5.1  has  a  few  advantages  over  Algorithm  4.1.  Algorithm  5.1 
does  not  incorporate  the  deflation  process  to  estimate  the  numerical  rank  as  often  as 
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Fig.  5.1.  Improved  Reduction  Steps  1-3  for  Downdating  the  ULVD 
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Fig.  5.2.  Improved  Reduction  Steps  4-6  for  Downdating  the  ULVD 
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Algorithm  4.1.  It  turns  out  that  after  a  downdate,  / the  ( k,k )  entry  of  1  often  gives 
a  very  good  approximation  to  Steps  3  and  5  offer  this  benefit  with  some  extra 

work,  namely,  2  (k  —  1)  Givens  rotations. 

Remark  5.3.  The  algorithm  can  be  also  used  for  updating  the  ULVD  just  by  replacing 
Step  4  with  the  Algorithm  3.8,  which  requires  only  12  flops. 


5.3  Rank  Detection 
5.3.1  Bounding  \\L  *|| 

In  this  section,  we  show  that  in  Step  8  we  can  often  determine  the  numerical  rank 
of  A  without  the  deflation  process.  We  start  with  stating  a  simple  lemma  that  bounds 
the  2-norm  of  a  block  matrix  in  terms  of  2-norms  of  its  blocks. 


Lemma  5.1  ([9,  Lemma  3.3]) 

.  Let  M  and  M  be  the  s  x  s  block  matrices , 

'Mn  Mu  - 

M  ) 

Is 

'  l|Mn|| 

iim12ii 

'•  IIM,/ 

M  = 

*  •  . 

,  M  = 

. 

KMS1  MS2  - 

M  j 

SS  / 

I|Ms2|| 

-  »MJIJ 

Then  \\M\\  <  ||M||. 

Proof.  See  the  proof  in  [9].  □ 

In  fact,  Lemma  5.1  holds  for  any  consistent  norms.  A  straightforward  application  of  this 
lemma  results  in  the  following  lemma. 
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Lemma  5.2.  Let 


v  l 

L.  K 


0  7/1 


(5.10) 


where  is  nonsingular  and  7  /  0.  Then, 


Hi-1,1  <  itfi  ’"‘lin-l 

1  0  7_ 1 


(5.11) 


Proof ’  It  is  easy  to  verify  that 


L-1  = 


-1,-1 
-7  Ln  w 


(5.12) 


Then  taking  norms  and  using  Lemma  5.1  yield  (5.11).  0 

The  following  lemma  shows  the  effect  of  Steps  3  and  5  of  Algorithm  5.1  on  l\!, 
(k,k)  entry  of  L^\  i  =  2,4. 


Lemma  5.3.  Let  L  be  defined  in  (5.10).  Suppose 


Lu  0 

L=  11  =  LQ 

~T  ~ 

\  w  7  / 


(5.13) 


where  Q  is  orthogonal.  Then, 


-1  -1,.  ,  llr- 1  ,,2x1/2 

=  7  (1  +  ll-^i  1  HI  )  • 


(5.14) 
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Proof.  From  (5.13)  we  see  that 


l~l  = 


_-l  ~T  j—\  ~-l 

-7  w  Ln  7 


where  the  expression  for  L  1  is  given  in  (5.12).  Thus  we  have 


rTz-1  = 


Comparing  both  sides  yields  (5.14).  0 

Combining  Lemma  5.2  and  5.3  gives  the  following  result. 


LEMMA  5.4.  Assume  the  hypothesis  of  Lemma  5.5.  Then 


ii^ii  ^(iiOi.r’.-r1), 


(5.15) 


where 


h  f1/  2  ,  2v  ,1  [7 T“  272  .  2  2 1 

i>(x,y,z)=  -{x  +y)  +  - y{x  +y  )  -4i  ^ 


(5.16) 


Proo/.  It  is  easy  to  show  that  for  any  2x2  upper  triangular  matrix, 


we  have 
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2fc\  1  /  2  -  *2  ,  2x  .  1  /,  2  .  ,2  .  2s2  .  2  2 

^(5)  =  -(a  +6  +c)+-y(o  +b  +  c  )  -4a  c  .  (5.17) 

Let  a  =  ||£~^||,  b  —  — j~ *||X~*u;||,  c  =  7  1 .  Then  by  Lemma  5.3,  we  see  that  +  c2  = 
7  ,  and  by  Lemma  5.2,  \\L  ||  <  a^S).  Hence,  taking  a  square  root  in  (5.17)  proves 

the  lemma.  □ 

The  following  lemma  states  necessary  conditions  to  ensure  a  correct  rank  estima¬ 
tion  for  the  ULVD. 

LEMMA  5.5.  Let  C  be  defined  in  (4-2),  and  let  r)  =  \\FL  1  |j.  Suppose 

||T_1||  ||G||  (1  +  t?)  <  tol,  77  <1.  (5.18) 

Therij 

ok{C)>tol>ck+l{C).  (5.19) 


Proof.  By  singular  value  interlacing  property  [63,  Theorem  7.3.9],  it  is  easy  to 
show  that 


>  <rmin(i)  =  l|£  'll  eM(C)  <  <7,(0)  =  ||<J||. 


(5.20) 


Thus, 

ak{C)  >  tol  >tol.j^>  11^11  >  <rk+1(C). 

This  proves  the  lemma.  0 

Thus,  the  decomposition  always  remains  rank-revealing  as  long  as  the  conditions 
in  (5.18)  are  satisfied.  Note  that  these  conditions  are  not  enforceable  unless  there  is  a 
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reasonable  gap  in  the  spectrum,  but  the  conditions  for  its  existence  are  weaker  than  the 
RRQR  conditions  of  Hong  and  Pan  [62]. 

A  cluster  of  singular  values  around  tol  would  cause  one  of  the  conditions  in  (5.18) 
to  be  violated,  and  the  rank  to  be  underestimated.  Since  we  keep  tracking  ||F||^,  and 
||G||^,  automatically  as  a  part  of  the  algorithm  as  illustrated  in  Section  5.3.3,  we  will 
always  know  when  this  happens. 

The  quantities  in  the  above  lemma  cannot  be  tracked  efficiently  using  our  algo¬ 
rithm.  However,  some  good  bounds  can  be  tracked.  That  leads  to  the  following  theorem 
that  uses  only  computed  quantities. 

Theorem  5.1.  Consider  the  downdating  procedure  in  Algorithm  5.1  and  the  related  up¬ 
dating  procedure  given  by  Remark  3. 5.  Let  C  be  the  matrix  before  updating  ( downdating ) 
and  partitioned  as  in  (4-3),  and  let  L  have  the  condition  estimate  k  ~  ||A  *||— *  such 
that 

°k{C)  >K>tol>  <rk+1(C). 

Let  C  be  the  matrix  after  updating  (downdating)  and  partitioned  according  to 


(5.22) 


where  ip(x,y,z)  is  the  function  defined  in  (5.16). 

If 
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(i)  k  >  tol 

(H)  ll-FII Fl*  =  V<  1 
(in)  ||G||f  <  tol(  1  +  7?)-1 

then 

as{C)>tol>as+1(C). 


If 


(iv)  K  >  tol 

M  '="V"f|lf  +  II^I||2-||/'esl|2  =  5  <  1 
M  yi|G|lJ.  +  l|f«sl|2  +  |fss|2  <  t°K  1  + 

then 

-  to1  -  as(C)- 


(5.23) 


(5.24) 


Proof.  Using  the  fact  that  for  any  matrix  A,  ||i4||  <  ||  A\\p,  it  is  simple  to  show  that 
the  conditions  (i)-(iii)  satisfy  the  hypothesis  of  Lemma  5.5.  Thus,  (5.23)  immediately 
follows. 

Let 

s-1  p-s+l  n—p 

l  ~ 

L  0  0 

C  =  F  G  0 
0  0  0 
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where 


4-1  1 


L  = 


L  0 
-T  J 


w 


ss 


5  —  1 


-  ?T 
,  w  —  L  e 


F  = 


/  -T\ 

W 


\  J  / 


•  F=l  *  L  )•  4  = 


G  = 


1  P-4 

V 

SS 

f~s  G 


1 

p-5 


Then,  it  is  simple  to  verify  that 


Ilf  Ilf  =  +  Nl2  -  ll/JI2 

iiS||J.  =  lieilp+li//+|fja 


Again,  the  conditions  (v)-(vi)  ensure  that  ||.F||^,  and  ||G||^,  satisfy  the  hypothesis  of 
Lemma  5.5.  Hence,  (5.24)  also  follows.  D 

The  conditions  (iii)  and  (vi)  make  certain  that  there  exists  a  reasonable  gap  in  the 
spectrum.  When  the  conditions  (i)— (iii)  are  satisfied,  updating  results  in  a  rank  increase. 
Similarly,  when  the  conditions  (iv)-(vi)  are  met,  downdating  results  in  a  rank  decrease. 
From  the  theorem  we  see  the  importance  of  keeping  \\F\\  as  small  as  possible.  Several 
ways  to  do  this  are  proposed  in  [41,  42,  104]. 


95 


The  theorem  still  holds  even  if  k  is  replaced  by  k  «  || L  !||f  in  (5.22).  We  can 
explicitly  compute  ||i  ||F  in  0(n  )  flops  as  described  in  the  foDowing  section.  In  some 

cases,  this  may  be  a  more  pessimistic  bound  than  the  one  in  Theorem  5.1. 

We  should  note  that  there  is  still  a  possibility  that  it  could  be  any  underestimate 
of  er^Z).  Thus,  if  the  conditions  of  this  theorem  are  not  satisfied,  then  the  condition  of 
L  will  still  have  to  be  estimated.  When  the  conditions  (iii)  or  (vi)  are  violated,  indicating 
that  there  exists  little  gap  between  them,  one  must  refactor  or  redefine  tol.  If  none  of 
these  work,  one  can  conclude  that  the  ULVD  may  not  be  suitable  for  tracking  subspaces 
for  the  problem  under  consideration. 

One  of  the  factors  that  might  affect  the  rank  estimation  using  this  scheme  is 
incorrect  rank  estimate  from  the  previous  update/downdate.  This  problem,  however, 
can  be  solved  by  computing  the  initial  factorization  using  the  SVD  for  an  accurate  rank 
estimation.  Although  the  SVD  is  more  expensive  to  compute  than  the  ULVD,  the  cost 
will  become  negligible  when  amortized  over  the  cost  at  the  subsequent  updating  and 
downdating  steps. 

Next  two  sections  will  show  how  to  track  efficiently  the  quantities,  ||L-1  ||F,  ||F||F, 
and  ||G'||f. 

5.3.2  Tracking  \\L~ 

We  begin  with  the  following  lemma. 

Lemma  5.6.  Assume  the  hypothesis  and  terminology  of  Lemma  5.3.  Then, 

l|i“1H2f.  =  l|i711|lf+3 

7 


(5.25) 
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Proof.  From  (5.12)  we  see  that 


II L 


-1, 


— 1  .,2 


=  Uu%  + 


i  + ||l 


Then,  (5.25)  follows  directly  from  Lemma  5.3.  □ 

Because  of  orthogonal  invariance  of  ||  •  ||^,  we  observe  that 

ii{£(3)r\  = 

||{X(5)}-1||f  =  ||{I<4)}-1||f 


where  we  used  the  terminology  of  Algorithm  5.1.  Thus,  the  steps  in  Algorithm  5.1  that 
we  need  to  consider  updating  ||X-1||^,  are  Steps  2,  4,  and  6.  The  following  lemmata  give 
the  formula  for  each  step. 

Lemma  5.7.  Assume  the  terminology  of  Algorithm  5.1.  Then 


IKi'Vl 


=  ii{i(3)r1nl+ 


kk  kk 


(5.26) 


Proof.  Since 


ii{i(3)r¥f.  = 
ii{£(4)r,i£  = 


ii{if,)}'1n2F+ 

« 

ii{i(,3)}'1n2f.+ 


1  +  ||{lP)}-,J3)||2 

<fi’>2 

i + ii{iS3)}-1»(3)n2 


(5.27) 


(5.28) 
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we  see  that 


IK-t'Wf. 


■(3K-1  (3)||2 


i  +  IKiS,'} 


i + ii{ii?r,«>(3)ii2 


€)2 


(5.29) 


By  Lemma  5.3  we  obtain 

i + ii{i(l)r1«’(3)ii2 

«l?>2 


<ig}2’ 


(5.30) 


Substituting  (5.30)  into  (5.29)  yields  the  result.  □ 

Since  is  not  available  at  Step  4,  and  ||{L^}  ^  ||^_,  =  ||{L^)}.  ^U^,,  we  update 
this  at  Step  5. 

Lemma  5.8.  Assume  the  terminology  of  Algorithm  5.1.  Then 


IK^VlIf 


IKi'Wf 


(1)\  1 1|2  llcH2 
)  Hr  {J(2))2 

(5.31) 

(5).-l,,2  ll<i||2 

>  h  {9(4),2 

(5.32) 

where  L^c  =  h,  and  {L^]^ d  =  ff' 
Proof.  Let 


{a(1)}t  *<? 


,  lW  = 


Lw  h 


o  9 


n  =  uTi^ 

(2)  3 


98 


Then  we  obtain 


ii{i(1)rVF 


iu£,2)r1iif 


. y’ 

<»ii> 

(»ii> 


(5.33) 


(5.34) 


—1  (2)  —1 

By  orthogonal  invariance,  we  see  that  ||{lA  '}  ||^  —  \\{LK  /  11^5  so  that 


ii{i(1)}'1ni+ 


i  +  IUi(1)}'T/,(1)ll2 


<*  T 


1  +  IU£<2)}-1A||2 


{5‘2)>2 


By  Lemma  5.3,  we  see  that 


(5.35) 


i + im(1)rT/i1)ii2 


(»s?)2' 


(5.36) 


Substituting  (5.36)  into  (5.35)  yields  (5.31).  The  proof  of  (5.32)  is  similar.  D 

Computing  || {Zr^2^}- 1 and  ||{L^6h_1|lF  requires  °(fc2)  floPs  for  forward  and 
backward  substitutions,  respectively,  and  ||{i^}  *||^,,  <9(1)  flops. 

Finally,  when  there  is  a  rank  change  in  Step  8,  we  use  the  relation, 


-  1  2  —12  l  +  ll'I>Ti111||2 

I*  +  {r  f  - 
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where 


\ 

0 

^ kk  / 


(5.37) 


2 

Thus,  it  takes  0(k  )  back  substitution. 

All  of  these  formulae  also  work  for  updating  as  well  as  downdating  for  obvious 

reasons.  Updating  and  downdating  procedures  differ  only  on  Step  4  when  updating  and 

(4) 

downdating  2x2  matrices.  Thus,  l is  computed  in  two  different  ways.  Therefore,  the 
formula  (5.26)  can  be  also  used  with  the  updating  procedure. 


5.3.3  Tracking  \\F\\F  and  ||<7||^, 

For  completeness  we  show  how  to  update  the  ||i?||^,  and  ||(7||^,  in  Table  5.1  al¬ 
though  it  was  partially  described  in  [13].  Here,  we  denote  =  [Cr^]e  ,t  =  3,4.  Note 
here  that  having  calculated  ||5^||  and  ||<7^||  to  compute  ||fr’^||^,,  computing  ||G^||f 
only  requires  0(1)  flops.  As  in  tracking  ||T-1||^,,  all  of  the  formulae  in  Table  5.1  also 
work  for  updating  as  well  as  downdating. 


Table  5.1.  Tracking  ||.F||  and  ||G||  for  the  Improved  ULVD  Procedure 

r  r 


Steps 

Updating  Formula 

Flops 

2 

IWinWWnBwnwM 

4 

wiaBffiMiiaBBmBMiaHa 

G(l) 

6 

BBSIItBIBilllBHIHilllBiiSlIM 

IBIHB1IHI11II1I 
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5.4  Numerical  Examples 

In  this  section,  we  present  a  few  examples  from  numerical  experiments.  These 

tests  were  performed  using  Matlab  on  a  SPARCstation  5  in  IEEE  Standard  double 

_ 16 

precision  with  machine  precision  «  10  .As  in  Chapter  4  the  algorithm  employs  the 

sliding  window  technique  from  signal  processing. 

At  each  step  of  the  sliding  window  method  with  the  window  size  an  mQ  x  n 

data  matrix  is  constructed  from  an  m  x  n  observation  matrix  A  by  adding  a  new  row  to 

the  data  matrix  in  the  previous  window  and  deleting  the  oldest  row  from  it.  In  step  j, 

the  row  +  j  of  the  observation  matrix  is  added  and  the  row  j  is  deleted,  giving  the 

data  matrix  A ..  The  ULVD  of  the  initial  window  matrix  An,  which  consists  of  the  first 
J  0 

rows  of  A,  can  be  obtained  by  computing  its  SVD. 

Then  Algorithm  5.1  takes  the  lower  triangular  matrix  (middle  part  of  the  decom¬ 
position),  the  orthogonal  matrix  (right  part)  as  initial  input  and  the  modifying  vector 

T 

r,  and  successively  modifies  these  matrices  at  every  window  step.  The  vector  z  =  V  r 
is  computed  at  the  beginning  of  each  window  step. 

We  tested  our  algorithms  in  the  context  of  the  total  least  squares  (TLS)  problems. 
See  Section  2.7  for  details.  We  use  the  TLS  solutions  from  the  Jacobi  SVD  as  reference 
in  checking  the  accuracy  of  the  solution  and  rank  estimates  of  our  algorithms. 

Fig.  5. 3-5. 5  show  the  rank  estimates  by  Algorithms  4.1  and  5.1.  The  horizontal 
axis  represents  the  window  steps  and  the  vertical  axis  the  numerical  rank  of  the  window 
matrix. 

The  distance  between  the  subspaces  is  given  in  the  next  plot  using  the  Definition 

a.  =  y.x.wt ,  W.  =  (W.  W.n) 

3  3  3  3  3  '  jl  3%' 


2.3.  Let 
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be  the  SVD  of  A .  computed  by  the  one-sided  Jacobi  method  at  step  j.  Let 


A .  =  V®  =  (V^  vj*)),  i  =  1,2 

3  3  3  3  3  3 1  3% 


be  the  ULVD  of  A.  computed  by  Algorithms  4.1  and  5.1,  respectively. 
Note  that  here  we  are  discarding  U\  ',  i  =  1,2.  Finally,  let 


c®  =  vHJ),  *  =  1,2 

3  3  3  1  3  J  3  x  3 1  J2  n 


by  the  SVD  of  C\  '  computed  by  the  one-sided  Jacobi  method.  Define  W .  by 


=  (  fi^W  i  =  1,2. 

3  K  j 1  j2  ’  3  3 


Define  the  angles  between  the  subspaces 


sin  ^  =||  {wWfvg  ||,  sin  if  =||  {W®)TV®  ||,  i  =  1,2. 


(5.38) 


The  angles  0^,  l  =  1,2  represent,  respectively,  error  between  the  true  noise  subspace  from 

the  Jacobi  SVD  and  the  approximate  one  from  tracking  the  ULVD,  the  approximation 

error  from  the  ULV  decomposition,  and  the  subspace  errors  from  ULV  subspace  tracking. 

We  plotted  log1Q(sin  0^)  in  dashed-dot,  log1Q(sin  0^)  in  solid,  and  log1Q(sin  0^) 

(2) 

in  dotted  line  on  the  vertical  axis  of  the  second  graph,  sin  0^  is  the  approximation 
errors  discussed  by  Fierro  and  Hansen  [43]. 
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Finally,  the  TLS  errors 


||l(SVD)  _  x(ULVD)|| 


i  =  1,2 


are  given  in  logarithm  in  the  last  plot.  Here,  and  i  =  1,2  are  the  TLS 

solutions  using  the  SVD  and  the  ULVD  with  Algorithms  4.1  and  5.1,  respectively.  On 
the  third  graph  of  each  figure  we  plotted  in  solid  and  in  dashed-dot. 

The  following  examples  were  also  used  in  [13]. 

Example  5.1.  A ,  a  110-by-6  random  matrix,  6,  a  110-by-l  random  vector.  Entries  of  A 

and  b  were  chosen  from  a  uniform  distribution  on  the  interval  (0, 1).  85  randomly  chosen 

—4 

rows  of  ( A ;  6)  were  multiplied  by  7  =  10  in  order  to  vary  the  rank  of  the  matrix,  and 
—2 

tol  =  10  .  The  window  size  p  used  was  12. 

—8  _ 6 

Example  5.2.  Same  as  Example  1  except  that  7  =  10  and  tol  =  10 

5 

Example  5.3.  Same  as  Example  1  except  that  the  matrix  had  an  outlier  of  size  10  at 
(18, 1)  position. 

The  first  plot  shows  that  both  algorithms  estimated  the  numerical  ranks  correctly 

throughout  the  sliding  window  steps  in  spite  of  frequent  rank  changes.  Thus  the  rank 

estimates  from  both  algorithms  and  those  by  the  Jacobi  SVD  were  identical.  As  expected, 

(i) 

the  errors  in  the  TLS  solution  are  almost  exactly  the  same  as  the  size  of  sin  0^  ,  i  =  1,2 
in  (5.38). 

The  second  plot  in  each  figure  shows  that  the  noise  subspace  error  is  very  small 
giving  accurate  TLS  solutions.  The  quantities  in  (5.38)  are  shown  to  be  essentially 
identical  indicating  that  the  subspace  errors  from  our  algorithm  are  from  the  rounding 
errors,  not  approximation  errors. 
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Moreover,  both  algorithms  perform  well  even  when  G  becomes  singular  (indicated 
by  for  Algorithm  4.1  and  for  Algorithm  5.1  in  the  first  plot).  We  tested  several 
other  examples,  and  these  results  were  typical. 

As  in  [13]  we  used  the  Corrected  Seminormal  Equation  (CSNE)  technique  when- 

o  n 

ever  the  downdating  is  not  possible,  namely,  ||a||  >  1  in  (4.5),  and  /3  +  /?  >  1  in  (5.6). 

1  2 

We  indicated  the  time  steps  where  the  CSNE  was  used  with  for  Algorithm  4.1  and 
’#’  for  Algorithm  5.1. 

For  Example  1,  we  plotted  the  norm  estimates  by  Algorithm  5.1:  ||T-1 1|  computed 
by  the  SVD  (in  solid),  R  of  Theorem  5.1  (in  solid  dot),  and  ||X-1||^,  computed  by  the 
Algorithm  5.1  as  described  in  section  5.3.2  (dotted),  again  on  a  log  scale,  relative  to  tol. 
This  verifies  the  bound  from  Theorem  5.1. 


Fig.  5.4.  Example  5.2 


Fig.  5.5.  Example  5.3 
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Chapter  6 


Modifying  the  Singular  Value  Decomposition 


6.1  Introduction 

We  discuss  methods  for  updating  and  downdating  the  SVD  and  partial  SVD  of 
A  6  7 ZmXn  of  the  form  (1.2)  and  (1.3).  Throughout  this  chapter  we  use  B  in  place  of 
M  in  (1.1)  to  denote  diagonal  form  or  partially  reduced  bidiagonal  form. 

Unlike  the  Jacobi-type  SVD  updating  procedures  [79,  80,  81],  we  transform  up¬ 
dating/downdating  problem  into  a  problem  of  finding  a  bidiagonal  matrix  B  and  an 
orthogonal  matrix  V  such  that 


BTB  ±  zzT  =  VBTBVT 
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We  may  also  use  the  MATLAB-like  shorthand 

B  =  bidiag(7(l:n),^(l:n-  1)) 

to  denote  the  above  bidiagonal  matrix.  As  in  modifying  the  ULVD  from  the  previ¬ 
ous  chapters,  our  approaches  to  downdating  the  decompositions  use  ideas  from  chasing 
algorithms  [1,  93, 115,  125]  and  from  the  downdating  algorithm  due  to  Saunders  [46,  85]. 
The  following  are  the  main  results  of  this  chapter: 

•  Procedures  for  updating  and  downdating  the  SVD  which  obtain  bidiagonal  forms 
such  that 


where  B  and  B  are  upper  bidiagonal,  and  /  =  k  +  1  for  updating  and  l  —  k 

1  Z 

for  downdating.  This  form  preserves  more  of  the  accuracy  of  the  small  singular 
values  and  is  not  achieved  by  standard  chasing  procedures.  We  can  then  use  one  of 
several  algorithms  to  find  the  singular  values  of  the  bidiagonal  matrix  B  to  relative 
accuracy  [11,  35,  40].  The  singular  vector  matrix  can  be  modified  by  a  procedure 
due  to  Gu  and  Eisenstat[54]  in  0(mn )  operations  (the  constant  on  mn  depends 
upon  machine  precision).  That  is  the  same  order  of  complexity  as  for  the  ULVD 
methods  with  similar  stability  properties. 

•  A  perturbation  theory  for  the  singular  subspaces  from  modified  matrices  and  block- 
wise  error  bounds  for  the  above  procedures. 
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The  condition  (6.2)  is  achieved  because  the  algorithms  for  both  updating  and 


downdating  problems  produce  an  orthogonal  matrix  V  that  has  the  form 


(6.3) 


There  is  never  a  rotation  of  the  first  k  columns  of  B  with  the  last  n  —  k. 

In  the  following  section  we  describe  the  secular  equation  approach  as  an  alterna¬ 
tive  to  chasing  algorithm.  We  show  that  this  approach  can  fail  to  separate  the  singular 
values  in  separate  blocks.  In  Section  6.3  we  give  some  basic  chasing  procedures  for 
modifying  the  SVD,  and  our  chasing  procedures  which  have  the  property  (6.2).  Section 
6.4  gives  the  perturbation  theory  and  discusses  error  analysis.  Section  6.5  gives  some 
computational  examples. 


6.2  Secular  Equation  Approach 

The  alternative  to  chasing  algorithms  for  modifying  the  SVD  is  that  of  finding 
the  zeroes  of  a  particular  spectral  function  [10,  25,  49,  53,  54,  67,  97], 

f(a)  =  1  +  a  Yh  '  2"'L"2  =  0  (6-4) 

i=l  ai  ~  d 

where  a  >  0  for  updating  and  a  <  0  for  downdating,  and  a  is  the  singular  value  of  the 
modified  matrix.  The  corresponding  singular  vector  is  given  by 


v  = 


,dTd  _2  T  s-1 

(BB  —  al)  z 
_2 _ n_ _ 

\\(BTB  -  o2lJ-'z\\ 


(6.5) 


no 


Here,  we  assume  that  B  is  diagonal.  That  approach,  as  yet,  does  not  allow  us  to  separate 
the  singular  values  into  separate  blocks  as  shown  in  the  following  example. 

Example  6.1.  Let 

1  0  'l 

B=  0  1 ;  10-10  ■ 

n  110~10; 

QR  decomposition  B  is  given  by  B  =  QR  where 


Q=  o 


-7.0711  •  10-1 

4.0825  •  10_1 

-5.7735  •  10-1 

0 

-8.1650 -10-1 

-5.7735  •  10_1 

-7.0711  •  10-1 

-4.0825  •  10-1 

5.7735  •  10-1 

-1.4142  -7.0711-10 


R  =  o 
,  o 


-1.2247-10 

0 


The  S VD  of  R  is  given  by  R  =  U  SV  where 


-1.0000  4.3301-10  0 

U  =  -4.3301- 10-21  -1.0000  0 

i  0  0  1.0000 


1.4142  0 

S=  0  1.2247  -10-10  >  v  = 

,00 


-5.0000-10 


5.0000-10  1 


Ill 


are  wrong! 


112 


6.3  Ordinary  Chasing  Algorithms 
6.3.1  Basic  Chasing  Routines 

In  this  section,  we  present  an  example  of  an  orthogonal  chasing  scheme  that 
produces  orthogonal  matrices  U,  V  6  TZnXn  such  that 

B  =  UTBV,  VTz  =  pev  p=\\z\\  (6.6) 

where  B  is  lower  bidiagonal.  For  the  4  x  4  case,  it  is  given  in  Fig.  6.1.  Here  r  and  x 
denote  possibly  large  elements  and  e  and  y  denote  small  elements.  See  Algorithm  3.4 
for  the  formal  description. 

In  theory,  this  could  be  used  to  produce  an  updated  or  downdated  bidiagonal 
matrix  very  easily.  We  have 


rp  rp  rp  _  *-*.rP  O  T1  —rP  — 

V  ( B  B±zz  )V  =  Bl  B±  p  CjCj  =  B  B 


If  BT  =  bidiag(7(l:n),^(l:n  -  1)),  then  BT  =  bidiag(7(l:  n),  n  -  1))  is  identical 
^rp  I  9  9 

to  B  except  that  7  =  —  p  .  This  is  illustrated  in  the  last  rotation  in  Fig.  6.1, 

denoted  by  a  pair  of  For  updating  it  is  simply  a  Givens  rotation.  It  should  be 
noted  for  downdating  that  the  assumption  |7j|  >  p  is  equivalent  to  the  assumption  that 
B  B  -  zz  is  positive  semi-definite. 

Such  a  procedure  shown  in  Fig.  6.1  does  not  preserve  the  separation  of  subspaces 
for  large  and  small  singular  values  as  accurately  as  we  would  like.  Large  elements  can 
get  chased  down  into  the  lower  part  of  the  bidiagonal  matrix  B  as  shown  in  the  following 
example. 
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11  11 


11  11 


Fig.  6.1.  Ordinary  Chasing  Procedure 
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Example  6.2.  Suppose  we  have 

B  =  diag(0.2071, 1.510  •  10-3, 8.081  •  10-4, 6.383  •  10-4, 5.184  •  10-7) 
z  =  (7.964  •  10-3, 8.012  •  10-3,  -9.102  •  10-3,  -2.821  •  10-3, 1.607  •  10-2)T 

After  forchase  is  applied  to 


we  obtain  B  =  bidiag(7(l:  5),  <^(1:4))  where 

7(1:  5)  =  (7.864  •  10~2,  -5.357  •  10-2,  -1.317  •  10-3,  -7.283  •  10-4,  -6.414  •  10_4)T 
^(1:4)  =  (-0.1852,4.134-  10-5, -4.030  •  10-4, 8.115  •  10-5)r 


Here  7(5)  >  tol,  so  that  the  smallest  singular  values  may  be  overestimated.  In  fact,  there 
should  be  no  rank  increase  as  we  will  see  in  Example  6.3. 

In  the  next  section,  we  show  that  forchase  and  backchase  procedures  described 
in  Section  3.2.1  can  be  combined  in  a  fashion  that  allow  us  to  update  or  downdate  the 
SVD  more  accurately. 

6.3.2  The  Updating  Algorithm 

Let 


B  =  diag(7(l:  n)),  7j  >  •  •  •  >  >  €  >  7]k+1  >  ' '  *  >  7„- 
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We  partition  B  into 


Bx  =  diag(7(l:  k)),  B2  =  diag(7(fc  +  1:  n)). 

Let  2  be  defined  in  (1.9).  Then  the  following  algorithm  updates  the  diagonal  matrix  into 
upper  bidiagonal  matrix. 

Algorithm  6.1  (Procedure  for  updating  diagonal  matrix).  Given  input  7(1:  n) 
that  contains  a.(A),i  —  l,...,n  and  the  update  vector  2,  this  procedure  produces  the 
updated  bidiagonal  matrix  B  =  bidiag(7(l:n),<^(l:n—  1)).  We  also  input  k  the  number 
of  singular  values  greater  than  tol. 

Step  1.  Compute  orthogonal  matrices  0  €  Kkxk  and  0  ,  V  € 

such  that 


0TlBlV  =  B^\  V'i  =  />,  =  11*11  (6.8) 

D2B2V2  =  B2)'  "2  =  M  I69) 

where  is  upper  bidiagonal,  and  B ^  lower  bidiagonal. 

1  Z 

Step  2.  Let  Q  =  J(l,3,0  )  and  Q  =  J(2,3,0  )  define  Givens  rotations  for  some 

X  X  Z  4 

0.,z  =1,2  such  that 
z 
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Note  that  after  Step  1  the  updating  problem  is  reduced  to  a  2  X  2  problem.  Step 
3  is  to  restore  the  matrix  to  upper  bidiagonal  form,  which  is  equivalent  to  a  step  of  qd 
procedure.  The  block  diagonal  form  of  V  in  (6.10)  is  highly  significant.  Any  modification 
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of  the  subspaces  associated  with  the  first  k  singular  values  and  the  last  n  —  k  singular 
values  will  be  computed  in  the  reduction  of  the  bidiagonal  matrix  B. 

Fig.  6. 2-6.3  show  the  reduction  steps  for  this  algorithm  with  n  =  7  and  k  =  4. 
Here,  a  pair  of  denotes  the  application  of  Givens  rotations  that  corresponds  to  Step 
2.  Unlike  the  ordinary  chasing  algorithm,  Algorithm  6.1  preserves  the  block  structure 
as  illustrated  in  the  following  example. 

Example  6.3.  When  Algorithm  6.1  is  applied  to  (6.7),  we  obtain 

7(1: 5)  =  (6.441  •  10" 4, 9.869  •  10-4, 3.974  •  10-3, 5.009  •  10-2, 3.818  •  10-7)T 
4>{l: 4)  =  (1.054  •  10-4, 7.309  •  10-4, 0.1862,  -4.444  •  10_8)T 

keeping  the  separation  of  the  large  and  small  blocks,  and  so  preserving  the  accuracy  of 
the  tiny  singular  values. 

6.3.3  The  Downdating  Algorithm 

An  immediate  difference  between  the  updating  and  downdating  procedures  is  that 
we  write  B  in  the  form 

H  =  diag(H1,B2,On_p)  (6.11) 

where 

B1  =  diag(7(l:fc))=  diag(<71,...,<Tifc), 

B2  =  diag(7(fc+l:p))  =  diag(<Tjk+1,...,<rp), 

”t>l>  Ck+r  V  l  =  -"  =  <,n  =  0' 


and 


118 


11  11  1111 


Fig.  6.2.  Bidiagonal  Reduction  Steps  for  Modifying  the  SVD 


119 


Fig.  6.3.  2x2  Updating/Downdating  Steps  and  a  qd  Step 
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thus  allowing  for  the  possibility  that  some  singular  values  are  exactly  zero.  As  in  down¬ 
dating  the  ULVD,  we  partition  the  vector  z  in  the  form, 


Here,  y  is  presumed  to  be  the  result  of  rounding  errors  and  is  ignored.  However,  if 
IIi/qII  >  M  *  ^(A),  then  we  should  not  downdate,  we  should  refactor. 

T 

Even  when  there  are  no  zero  singular  values,  and  even  though  z  —  V  r  and  r  is 
T  T  . 

the  first  row  of  A  in  (1.1),  B  B  —  xx  is  positive  semi- definite.  However,  occasionally, 

T  T 

even  that  is  not  the  case.  The  usual  way  to  test  if  B  B  —  xx  is  positive  semi-definite 
is  to  solve 

B^s  =  x.  (6.12) 

If  \\s\\  >  1,  then  we  cannot  downdate  B ^  by  x.  One  possible  remedy  is  to  try  to  obtain 

T  T 

a  better  value  for  x  =  V  A  e^.  That  can  be  done  using  the  corrected  semi-normal 
equations  (CSNE)  [18,  20]  as  we  have  used  for  modifying  the  ULVD  in  Chapter  4. 

If  ||*||  >  T}  where  t?  <  1,  then  solve 

BlC  =  s.  (6.13) 

It  should  be  noted  that  c  solves  the  least  squares  problem 


min  \\AVc-e  || 
c€Tl 
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and  its  value  can  be  improved  by  the  iterative  improvement  steps 


r  =  el  -  Avic 

(6.14) 

6x  =  ATr 

(6.15) 

B^6s  =  6x 

(6.16) 

5  +  6s,  X  X  +  Sx 

(6.17) 

At  this  point,  if  ||s||  >  1,  we  signal  that  downdating  is  not  possible,  and  thus  other  options 
should  be  considered,  such  as  refactoring  or  choosing  a  higher  threshold  e.  Otherwise, 
the  algorithm  proceeds  in  a  similar  manner  to  Algorithm  6.1. 

We  now  present  the  downdating  algorithm. 

Algorithm  6.2  (Procedure  for  downdating  diagonal  matrix).  Given  the  input 
7(1:  n )  that  contains  a. (A), i  =  1, . .  .,n  and  the  update  vector  z  of  the  form  (4.3),  this 
procedure  produces  the  downdated  bidiagonal  matrix  B  —  bidiag(7(l:»),<^(l:n—  1)). 
We  also  input  k  the  number  of  singular  values  greater  than  tol.  yQ  is  ignored  unless 

\\y0\\  >  M* 

Step  1.  If  ||s||  >  tj  where  s  is  defined  in  (6.12)  and  r)  <  1,  then  solve  (6.13)-(6.16) 
for  6s  and  6x.  Update  s  and  a:  as  in  (6.17).  If  ||s||  >  1  or  ||i/q||  >  tol  a  100  *  p, 
then  quit  and  exit;  otherwise,  do  Steps  2-4. 

Step  2.  Same  as  Step  1  of  Algorithm  6.1. 

Step  3.  Compute 


(6.18) 


a  = 


(6.19) 
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We  can  then  find  two  3x3  Givens  rotations  Q .  =  J(l,i  +  1, 9.),  i  =  1, 2  such  that 


^1^2  ai 


In  that  case  we  modify  the  Ar-th  and  ( k  +  l)-st  rows  of  the  matrix 


42)  42) 

o  T« 


7^  0 


=  «1«2  0  iff,  ' 


0  0 


where  p ^  =  7^^  V  1  —  <?•  This  can  be  done  by  using  Algorithm  3.9, 


^2)’42)^Sll  =  down22(7i1),°,7^i,ai,a2). 


Thus  if  we  let  =  J(l,k  +  2, 6^  J(  1,  fc  +  1,  then  we  have 


n(2)  (2)  T 

B\  4  Vi 
0 


V*  Vi 


'if*  0 

=  ^  0  if 

,  0  0 


Step  4.  Same  as  Step  3  of  Algorithm  6.1. 


We  note  that  in  (6.18)  laj  =  ||j x||  =  ||s||.  Thus  B ^  can  be  downdated  by  x  if  and 

only  if  |a  |  <  1.  For  Algorithm  6.2,  we  assume  that  is  the  case.  If  a  =  0,  but  a  /  0, 
1  ^ 
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V*  ('’2  +  ^1 

5i  hv* 


(6.23) 


Thus  from  (6.20),  ||£r||2  =  6p 2  +  £p2.  Using  (6.23)  and  noting  that  Ue ^  =  Cj,  we  have 


(T  F  = 

5  i4  +  «i4 


(6.24) 


Thus  comparing  (6.24)  with  (6.21),  we  have  6  A. 


foUows.  □ 


(o  aT). 


The  result  immediately 


We  note  that  Proposition  6.1  is  merely  a  consistency  property.  What  it  says  is 
that  approximation  used  in  Step  3  of  Algorithm  6.2  does  not  increase  the  error  over  that 
caused  by  assuming  that  Ue^  =  e^. 

6.3.4  Extensions  to  Partially  Reduced  Bidiagonal  Forms 

Algorithm  6.1  and  6.2  can  be  easily  extended  to  the  case  where  either  B  or  B ' 

is  bidiagonal  as  long  as  they  are  decoupled.  We  need  only  modify  Step  1  of  Algorithm 

6.1  and  Step  2  of  Algorithm  6.2.  Van  Huffel  and  Park  [115]  describe  chasing  algorithms 

kxk 

that  given  B ^  upper  bidiagonal,  produce  orthogonal  matrices  such  that 


where  B ^  is  upper  bidiagonal. 

Fig.  6.4  shows  how  such  an  algorithm  would  work  on  a  5  x  5  example.  Thus, 


using  algorithms  such  as  the  zero-shift  QR  [35]  or  the  qd  algorithm  [40],  it  is  possible  to 
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Fig.  6.4.  Reduction  Steps  for  the  Partially  Reduced  Bidiagonal  Form 
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find  the  singular  values  that  are  below  a  certain  threshold,  and  thus  obtain  a  partially 
bidiagonal  matrix  of  the  form  (1.3). 

6.4  Error  Analysis 

6.4.1  Error  Bounds  for  Blockwise  Algorithms 

We  now  present  error  bounds  for  the  process  of  one  update  or  downdate  using  the 
procedures  in  Sections  6.3.2  and  6.3.3.  All  of  the  matrices  below  are  computed  except 
those  with  6  in  front  of  them. 

The  following  two  propositions  are  proven  in  the  Appendix  of  [14]. 

Proposition  6.2.  Algorithm  6.1  produces  an  updated  matrix  B  such  that  for  some  or¬ 
thogonal  matrix  U,  and  V  we  have 


T  \ 

(  „  \ 

-T 
U 1 

z 

K  B  ) 

V  = 

0 

^  B  +  SB  +  SBq  j 

where 

6B  =  di*g(6Bv6B2) 

sBo  =  lVtet  +  %n  vf+i +  ,T*+Ie*+Iet+1 

1*5,11  <(*/,(»)  l|5,ll  +  <V) 

||*5jll  <  c/2(")  P2II  +  OOO  =  (*/2W  °m(B)  +  Oin1) 
1*7^1  <  m/3(")  1 7,4  +  0((i2).  j  =  k,  k  +  1 

l**t+1l  <*•/„(»)  l^+1l  +  0("2) 

where  /,(n)  =  0(r?),  i  =  1,2,  and  /,(n)  =  0(n),  i  =  3,4. 
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From  [11,  35]  this  relative  change  in  the  entries  of  the  bidiagonal  matrix  makes 
only  relative  changes  in  the  singular  values.  Thus  this  update  procedure  is  very  stable. 

PROPOSITION  6.3.  Algorithm  6.2  produces  a  downdated  matrix  B  such  that  for  some 
orthogonal  matrices  U  and  V  we  have 


\ 

,  p  <r_  \ 

B  +  6B 

z  +  szyv 

v  = 

K  0  ) 

l  *  i 

where 


SB  =  diag  (SBvSB2) 

n«yi  < /*/5(»)  11^11+ <V) 
ll**2ll  <  M/e(»)  \\S2\\  +  <V)  =  m/6(»)  °k+l{B)  +  <V) 

||<5z(l:  k) ||  <  (if7{n)  ||z(l:  k)\\  +  0(^i2) 

||^z(fc+  l:p)||  <  nfg(n)  \\z(k+  l:p)||  +  0(fi2) 

where  f^n)  =  0(n2),  i  =  5,6,  and  f^n)  =  0(n),  i  =  7,8. 

These  results  are  as  good  as  can  be  expected  for  any  such  procedure.  As  we 
state  in  the  next  section,  we  can  expect  sharp  separation  between  singular  subspaces 
associated  with  large  and  small  singular  values. 

6.4.2  Perturbation  Bounds  for  Invariant  Subspaces 

We  consider  in  this  section  the  effects  of  the  bounds  in  Propositions  6.2  and  6.3 
in  the  error  in  certain  invariant  subspaces  of  B  resulting  from  Algorithms  6.1  and  6.2. 
Two  perturbation  results  show  that  we  expect  that  the  subspaces  for  large  and  small 
singular  values  will  be  very  accurately  computed. 
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The  componentwise  backward  error  6BQ  in  Proposition  6.2  has  a  very  small  effect 
on  the  error  in  the  subspaces.  The  following  result  is  given  for  completeness. 


Proposition  6.4  ([39,  Lemma  4.5]).  Let 


B  =  bidiag(71,...,7n;^1,...,<^n_1) 


(6.25) 


and  let 


B  =  bidiag(al7l ,  •  •  • ,  «n7n;  «n+1^ .  •  •  • » a2n-\K-\^  ^ 


7=  I]  maxima.  }  -  1. 
i=  1  ~  J 


Let  be  the  right  singular  vectors  of  B  and  let  5^ , . .  •  be  the  right  singular 

vectors  of  B.  Let  <7  , . . . ,  <r  be  the  singular  values  of  B  and  define 


f  K  -  °i\  \ 

p.  =  min<2,min — : — : — >,  i=l,2, 

1  1  i#»  K-l 


Let  Z.  =  (w.,...,w.  ,,w.  u>  ),  that  is,  the  right  singular  vector  matrix  of  B 

i  v  1’  ’  1— 1’  i+l  n ' 


with  its  i-th  column  deleted.  If  p ^  >  r),  then 


11  V  p,  -7?  2J 


(6.27) 


Thus  the  effect  of  the  relative  errors  6BQ  on  the  updated  matrix  B  is  minimal 
and  has  little  effect  on  the  singular  subspaces. 
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Proposition  6.3  implies  that 


I|w,II<*bI|b,II  +  o((.2) 

(6.28) 

II**2II<  V  +  o(m2) 

(6.29) 

l|fa(l:*)ll£*,Wl:*)ll  +  0(»<J) 

(6.30) 

[|J+(fc  +  1:  p)||  <  b}\z(k  +  1:  ?.)||  +  Oi,i2). 

(6.31) 

Here,  6  <  f  (n)p  and  S  <  f  ( n)p  where  fn{n)  =  0(r, t2)  and  /  ( n )  =  0(n). 

±J  LJ  Z  Z  ±J  Z 

Proposition  6.5.  Let  B  and  B  +  SB  be  diagonal  matrices  such  that 


k  n-k 


k  n—k 


-  -  t\ 

bi  <**+1***1 

0  n—k 


SB  = 


k 

n—k 


where 


II**, II  <  *BP,II  +  °(A  II«*2II  <  V  +  0(1?), 


Let  a,  >  •  •  •  >  <j,  >  e  >  a,  ,>•••>  a  .  Let  W,  W  G  TZnxn  be  the  matrices  of  right 
1  —  ~  k  k+1  —  —  n 

singular  vectors  of  B  and  B  +  SB,  respectively.  If 


W  =  (W1  w2),  W  =  (Wl  w2), 


(6.32) 
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where  W  ,  W  €  ft  ,  and  W  , 


nxk,  and  W  ,  Wn  €  ftnx("  k\ 


\\^W2\\F  <  2^11^11  ||5j  l\\^±^l  +  0(\\6B\\2). 

k  A?4~l 


(6.33) 


Proof.  Let  uL,  i  =  1,2, . . k  be  the  i-th  right  singular  vector  of  B  +  SB  and  let 
w j  =  k  +  1, . . . ,  n  be  the  j- th  singular  vector  of  B.  Then  from  standard  perturbation 

-T  - 

bounds  on  the  eigenvectors  of  B  B  we  have 


\vk  Bw .  +  wTb^SBw.\ 

IffiVl  =  — - i — 5 - —  +  0(||$ij||2) 

J  a.  —  a . 

*  J 


\\6Bw.\\o.  +  \\SBw  \\a 


2  2 

a .  -  a . 

x  j 


+ mm 


(6.34) 


(6.35) 


We  now  bound  ||^5u;.||  for  i  =  1,2, ...,n.  First,  let 


w.  = 

*  ,  J 2) 


(6.36) 


Then  we  have 


*i  V i 
0 


(!)  \  (1)  \ 

w:  1  y)  ' 

1  —  o.  1  =  a.y. 

(2)  *  (2)  *  1 

Wi  \  yi 


where  y.  is  the  corresponding  left  singular  vector.  Thus  we  have 


(1)  n-\,  (1)  7  (2L 

wi  =  -®i  M  i”i  >• 
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Therefore, 


ihS1)n=pr,ii(<’i+ 1^1^11^11^+^ 


(6.37) 


Now  we  can  say  that 


||«».||2  <  ||WI»<1)||2  +  ||«52i»f  >||2 


which  leads  to 


<  «Bp,||  p,  'M^  +  ef  +  llc2. 


(6.38) 


Equation  (6.38)  leads  to 


ll*fi«{ll£V<’i  +  ^OP1ll  Ill'll- 


(6.39) 


Combining  (6.35)  and  (6.39)  yields 


T  __i  (<Ti  +  v/2  0  *•  +  (*,+  V2  e)  <ri 

I«f*-I  <  ^ll^ll  Pj  II— - i — r - 1  +  0(ll55ll ) 

cr.  —  <7  . 
z  j 


which  is  bounded  by 


7 ■»  _ i  (cr.  +  \/2  e)  _  9 

i«; *.i < y  ^11  HBj ■  +  odi^n2). 

*  i 


(6.40) 


Thus  for  all  i  =  1, 2, . . . ,  fc  and  j  =  fc  +  1, . . . ,  n,  we  have 


T  -  1  (^i.  i  i  -  2 

|ef»,l  <  «BIP,II  IP,  ll-^-T - +  o(||H)||2). 

3  ak~ak+ 1 


(6.41) 
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which  is  now  independent  of  i  and  j.  The  use  of  the  Frobenius  norm  on  =  (w^  ,...,w^) 
and  W  =  (w,  , . . w  )  yields  (6.33).  □ 

Z  K  “j- 1  7Z 

Proposition  6.5  implies  that  the  updating  algorithm  will  always  yield  accurate 
subspaces  for  the  first  k  and  last  n  -  k  singular  Values.  For  Algorithm  6.2  we  must  also 
bound  the  effect  of  8z  which  is  qualitatively  slightly  different. 

~T  - 

Proposition  6.6.  Let  B  and  W  be  an  in  Proposition  6.5  and  let  it  satisfy  B  B  — 
V iF BV^  —  zF  where  V  £  lZnxn  is  the  orthogonal  matrix  from  Algorithm  6.2.  Let  B 
satisfy  pF B  —  V B^ bF  —  (z  +  6z)(z  +  6z)^  and  let  8z  have  the  form  (6.30)-(6.31). 
Let  z  =  {F  FF >  x  £  ,  y  £  TZn~ ^  and  assume  that  ||i/||  <  e.  Let  W  £  TZnxn  be  the 
matrix  of  right  singular  vectors  of  B  and  define  and  from  (6.32).  Then 

II w(w2\\F  <  «^*(«-*).uli||B1||  llfi-'ll  x 

<r,  +  3c  4€2  o 

-  *  +  ~2~~ . —  +  0(\\6z\\2)  (6.42) 

k  ak+ 1 

where 

u;  =  \\B~Tx\\,  k  =  Il-Bjll  \\B~l\\  (6.43) 


Proof.  Let  w.,  i  =  1, 2, . . . ,  k  be  the  i- th  right  singular  vector  of  B  +  SB  and  let 
w . ,  j  =  k  +  1, . . . ,  n  be  the  j-th  singular  vector  of  B.  Then  from  standard  perturbation 

-T  - 

bounds  on  the  eigenvectors  of  B  B  we  have 


~T  , 
w .  w  A 

i  3 


+  o(IMI2) 


(6.44) 
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\SzTw  \  Iz^TO  |  +  \SzTw  \  \zT W  I 

<  - i - j— 1 - L  +  0(INn-  (6.45) 

a.  -  a . 

*  3 


If  we  use  the  partitioning  of  w.  in  (6,36),  then  we  have 


\zTw.\  =  +  \yTw{2)\  =  I xTB  lB  w\l) |  +  |/^2)| 

2  2  2  12  2 


which  means  that 


\z^w.\  <  || B  ^x||  +  €  <  a?  o\  +  €,  i  =  1, . .  .,n.  (6.46) 


We  also  have  that 


|far®.|  <  >1  +  l«/«f 'l  <  «,H  ll>»,.1)ll  +  Ml  ||«f  >||. 


Using  the  facts  that  ||x||  <  and  ||i/||  <  e  and  from  (6.37),  we  have 


|foT-»,l<<,(P1ll  llViK+f)+f 


which  we  simplify  to 

\6zTw.\  <  6JIBJI  \\B;l\\(<T.  +  2e).  (6.47) 


Combining  (6.45)  with  (6.46)  and  (6.47)  yields 


T  _i  K-  +  2<0K-  + e)  +  (*•  +  2e)(fri  + €) 

\wjw  \  <  6»\\B\\\\B  /H-i - 1 - 2 - T - 

1  J  *  1  1  0,-0. 

*  3 

+  o(IIMI2) 


(6.48) 
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< 


Ill'll 


( 7  .  -j-  36 
3 _ 


a.  —  a . 
*  3 


+  o(||MI2). 


(6.49) 


We  note  that  this  for  i  =  1, 2, . . . , k  and  j  =  k  +  l,...,n,  thus  we  can  bound  (6.49)  by 


\wjwj\  <  Su  HBJ  || Bl  !|| 


_VtfL  +  _iL_ 

ak  ~  ak+ 1  o\  -  v2k+1 


+  0(||^||2). 


(6.50) 


— j'  ~ 

The  bound  (6.42)  is  obtained  by  computing  the  Frobenius  norm  of  W '  where  = 

(®1»  •  *  *’«*)  and  W2  =  (wk+V  •  D 

6.5  Numerical  Examples 

In  this  section,  we  present  a  few  examples  from  numerical  experiments.  These 

tests  were  performed  using  Matlab  on  a  SPARCstation  5  in  IEEE  Standard  double 

precision  with  machine  precision  «  10  .As  in  Chapter  4  the  algorithm  employs  the 

sliding  window  technique  from  signal  processing. 

At  each  step  of  the  sliding  window  method  with  the  window  size  mQ,  an  mQ  x  n 

data  matrix  is  constructed  from  an  m  x  n  observation  matrix  A  by  adding  a  new  row  to 

the  data  matrix  in  the  previous  window  and  deleting  the  oldest  row  from  it.  In  step  j, 

the  row  m  +  j  of  the  observation  matrix 

Then  Algorithms  6.1  and  6.2  take  the  diagonal  matrix  and  the  orthogonal  matrix 

(right  part)  as  initial  input  and  the  modifying  vector  r,  and  successively  modifies  these 

T 

matrices  at  every  window  step.  The  vector  z  =  V  r  is  computed  at  the  beginning  of 
each  window  step. 
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The  one-sided  Jacobi  method  [36]  was  used  to  compute  the  SVD  of  the  initial 
window  matrix  which  consists  of  the  first  m  rows  of  A,  and  to  compare  with  our 
algorithms  for  rank  estimation  and  the  accuracy  of  the  subspaces. 

We  tested  our  algorithms  in  the  context  of  the  total  least  squares  (TLS)  problems. 
See  Section  2.7  for  details.  We  used  the  TLS  solutions  via  the  Jacobi  SVD  as  reference 
in  checking  the  accuracy  of  the  solution  and  rank  estimations  of  our  algorithms. 

In  Fig.  6.5-6. 7,  the  rank  estimated  by  our  algorithms  (solid  line)  and  the  true 
rank  (dotted  line  but  not  visible  in  the  plot)  are  given  in  the  first  plot.  The  horizontal 
axis  represents  the  window  steps  and  the  vertical  axis  the  numerical  rank  of  the  window 
matrix. 

Let  and  be  the  right  singular  vector  matrices  computed  by  the  Jacobi 
method  and  Algorithms  6.1  and  6.2,  respectively.  Then,  using  the  Definition  2.3, 

6.  =  \\W^TV^\\,  j  =  1,2, ... 

where  =  (W^  W^)  and  =  (V^  V^).  We  plot  log1Q(sin(^.))  in  the 

second  plot  of  each  figure. 

Finally,  the  TLS  errors 


are  given  in  logarithm  in  the  last  plot.  Here,  x ^  and  x  .  are  the  TLS  solutions  using  the 
Jacobi  method  and  our  algorithms,  respectively. 

Example  6.4.  A,  a  100-by-5  random  matrix,  6,  a  100-by-l  random  vector.  Entries  of  A 
and  b  were  chosen  from  a  uniform  distribution  on  the  interval  (0, 1).  75  randomly  chosen 
rows  of  [ A ;  b ]  were  multiplied  by  7  =  10-4  in  order  to  vary  the  rank  of  the  matrix,  and 
tol  =  10  .  The  window  size  p  used  was  10. 
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The  first  plot  shows  that  our  algorithms  estimated  the  numerical  ranks  correctly 
throughout  the  sliding  window  steps  in  spite  of  frequent  rank  changes.  The  errors  in  tiny 
singular  values  were  relatively  large,  and  our  algorithms  almost  always  overestimated 
small  singular  values.  However,  they  were  close  enough  for  the  correct  rank  estimation. 

The  second  plot  in  each  figure  shows  that  the  noise  subspace  error  is  very  small 
giving  accurate  TLS  solutions. 

—9  —7 

Example  6.5.  Same  as  Example  6.4  except  that  7  =  10  and  tol  =  10 

5 

Example  6.6.  Same  as  Example  6.4  except  that  the  matrix  had  an  outlier  of  size  10 
at  (15, 1)  position. 

Both  TLS  solution  errors  and  the  noise  subspace  errors  show  that  our  algorithms 
give  very  accurate  approximation  to  the  subspaces  under  consideration.  Moreover,  the 
algorithm  performs  well  even  when  some  of  tiny  singular  values  become  almost  zero 
(indicated  by  in  the  first  plot).  We  tested  several  other  examples,  and  these  results 
were  typical. 

Since  our  downdating  procedures  use  LINPACK  downdating  algorithm,  it  is  not 
difficult  to  generate  the  cases  where  the  algorithm  breaks  down  when  ||a||  >  1,  for 
instance,  when  deleting  a  large  row  relative  to  7  (see  Fig.  6.6)  or  a  row  that  contains 
outliers  (see  Fig.  6.7).  We  used  the  CSNE  approach  in  (6.14)-(6.17),  and  indicated  these 
steps  by  ’+’  in  the  first  plot. 

The  CSNE  approach  was  used  in  all  three  examples  and  most  extensively  in 
Example  6.6  when  downdating  a  row  with  an  outlier.  However,  the  performance  of  our 
algorithm  was  less  satisfactory  for  the  larger  outlier. 
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Chapter  7 


Parallel  Implementation 


7.1  Introduction 

In  this  chapter  we  describe  a  fully  parallel  bidiagonal  reduction  procedure  for 
modifying  the  SVD.  In  Chapter  6  we  showed  that  blockwise  algorithms  produced  more 
accurate  subspaces  than  the  ordinary  one-way  chasing  algorithms  which  ignore  the  block 
structure  of  the  diagonal  matrix.  A  VLSI  implementation  of  the  similar  chasing  schemes 
for  the  bidiagonal  reduction  for  updating  was  also  described  in  [1,  117],  but  without 
considering  the  large  and  small  structure  of  the  matrix.  In  this  section  we  implement 
our  algorithm  on  a  distributed-memory  MIMD  multiprocessor.  Two  storage  schemes  are 
considered:  cyclic  storage  scheme  and  consecutive  storage  scheme.  We  will  show  that  the 
consecutive  storage  scheme  implements  the  bidiagonal  reduction  much  more  efficiently. 

The  main  idea  behind  the  Algorithms  6.1  and  6.2  is  to  reduce  the  entries  of  the 
vectors  x  and  y  in  opposite  order  and  to  chase  the  bulge  in  opposite  direction,  upper-left 
corner  for  the  large  block  and  lower-right  corner  for  the  small  block.  This  is  based  on 
the  two-way  chasing  scheme  [125],  which  was  also  used  in  [116]  with  k  =  n/2  in  the 
context  of  updating.  The  algorithm  simply  reduces  the  large  and  small  blocks  to  almost 
bidiagonal  form  (see  the  12-th  matrix  in  Fig.  6.2)  by  ordinary  chasing  scheme,  and  uses 
2x2  updating  and  downdating  algorithms  to  eliminate  x ^  and  j/^ ,  followed  by  one  step 
of  the  qd  procedure  on  the  small  block  to  reduce  it  to  the  upper  bidiagonal  matrix.  The 
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entire  reduction  steps  for  Algorithm  6.1  requires 

Step  1  Step  3 

k(k  —  1)  +  (n  —  k)(n  -  k  —  1)  +  (n  -  k  —  1) 

=  n  -  2 nk  +  2 k2  -  (k  +  1)  (7.1) 

plane  rotations. 

The  algorithm  allows  simultaneous  bidiagonal  reductions  on  both  large  and  small 
blocks,  B  and  B  (Step  1  in  Algorithm  6.1)  since  they  do  not  share  any  data  throughout 

1  JU 

the  reduction  steps.  Following  similar  notations  used  in  [1],  we  denote  F.  .  as  Givens 

plane  rotations  operating  on  rows  i  and  j  (left  rotations),  and  G .  .  as  those  operating 

i)3 

on  columns  i  and  j  (right  rotations).  Then  from  the  dependency  graph  of  this  algorithm 
depicted  in  Fig.  7.1,  we  see  that  G  and  G  ,  the  first  rotations  for  each  block,  can 

i  ljTl 

be  executed  in  parallel,  and  the  sequence  of  the  rotations  that  follows  are  also  carried 
out  in  parallel.  Hence,  the  whole  reduction  only  takes 

k—Z 

2  +  3  +  —  +  3 +2(fc  -  1)  =  5k  -  9  if  k  >  (7.2) 

k- 3 

2  +  3  +  •  -  +  3  +2(n  —  k  —  1)  +  1  =  5n  —  5k  —  8  if  k  <  L^rJ  (7.3) 

time  steps.  Obviously,  the  algorithm  achieves  an  optimal  performance  when  k  ~  n/2. 

7.2  Overview  of  Connection  Machine 

A  Connection  Machine  (CM-5)  system  can  have  up  to  16K  physical  processors 
or  processing  nodes  (PNs).  The  CM-5  has  two  interprocessor  communications  networks: 
data  network  and  control  network.  The  control  network  is  used  for  global  operations 


s; 
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such  as  synchronization  and  broadcasting.  The  data  network,  which  uses,  so-called,  J- 
ary  fat  tree  [72],  supports  operations  for  data  transfers  from  a  single  source  to  a  single 
destination. 

The  CM-5  supports  both  SIMD  (Single  Instruction  Multiple  Data)  and  MIMD 
(Multiple  Instruction  Multiple  Data)  programming  models  [66].  In  the  SIMD  model,  the 
data  parallel  programming  associates  one  PN  with  each  element  of  a  data  set.  All  PNs 
execute  identical  operations,  each  operating  on  data  stored  in  its  local  memory,  accessing 
data  stored  in  the  local  memory  of  other  PNs,  or  receiving  data  from  the  host  computer. 

In  the  MIMD  programming  model,  each  node  has  its  own  copy  of  the  same  pro¬ 
gram  called  node  program,  and  executes  the  program  asynchronously.  The  communi¬ 
cation  between  the  nodes  is  usually  done  by  utilizing  a  set  of  efficient  communication 
routines  contained  in  the  CM  message-passing  library,  CMMD.  For  our  implementation 
we  chose  the  MIMD  model  because  the  rotations  at  each  time  step  are  different  in  terms 
of  their  types  (left  or  right)  and  the  data  required. 

In  CM-5  a  packet  of  size  20  bytes  is  used  for  nodal  communications.  First  four 
bytes  are  used  for  control  purposes  and  the  rest  of  16  bytes  contain  the  data.  If  a 
packet  is  full,  that  is,  if  it  contains  16  bytes  of  user  data,  the  overhead  of  processing  it  is 
smaller  than  the  message  of  different  sizes.  Therefore,  the  communication  overhead  will 
be  smaller  if  a  user  made  the  message  size  a  multiple  of  16  bytes  [90]. 

Moreover,  a  cluster  is  composed  of  four  processing  nodes,  and  the  nodes  with  the 
same  cluster  share  a  common  switching  node  capable  of  four  times  the  bandwidth  of  the 
node  at  the  leaf  level.  A  Similar  statement  is  true  for  the  nodes  as  progressing  toward 
the  root.  Each  node  must  go  through  at  least  one  switching  node  to  communicate  with 
the  other  node.  To  communicate  with  the  node  in  a  different  cluster,  the  communication 
path  will  be  longer.  Hence,  it  takes  longer  to  transfer  the  data  to  the  node  within  the 
same  cluster  than  to  the  node  outside  the  cluster. 
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We  assume  the  CM-5  consists  of  p  processing  nodes,  and  denote  them  as  Node(O), 
. . . ,  NoDE(p—  1).  Here,  we  assume  that  p  is  a  power  of  two.  We  choose  one  of  the  CMMD 
programming  styles,  Host/Node  model,  where  the  host  processor  allocates  the  data  to  the 
nodes,  and  collects  the  results  for  the  analysis.  The  host  program  calls  the  node  program 
residing  in  each  node  for  the  various  tasks,  and  each  node  has  identical  node  program. 
Once  the  node  program  is  loaded  in  each  node,  it  can  be  executed  asynchronously.  Both 
host  and  node  programs  for  our  implementation  are  written  in  FORTRAN. 

7.3  Implementation  Details 

In  this  section  we  give  a  detailed  description  of  parallel  implementation  of  the 
Algorithm  6.1.  First,  we  need  the  following  definition. 

DEFINITION  7.1.  A  pair  of  left  rotations  F.  . , ,  and  F.  .  .or  riqht  rotations  G .  .  . 

y  J  J  M+l  J,J+l  *  i,i+l 

and  G .  ._|_i  is  said  to  be  adjacent  if  | i  —  j  \  =  2. 

We  also  use  Lred(i)  to  denote  the  sequence  of  plane  rotations  for  eliminating  x. 
and  restoring  the  resulting  matrix  into  the  upper  bidiagonal  matrix  (large  block),  and 
Sred(i)  to  denote  those  for  eliminating  and  restoring  the  matrix  into  the 

lower  bidiagonal  form  (small  block),  that  is, 


LRED(t)  =  {G 
SRED(i)  =  {G 


F  F  G 

t,t+l’  i.i+1’  ’  1,2’ 

F 

n—i,n—i+ 1’  n—i,n—i+ 1’ 

i  =  1, —  fc  — 


12>,  i=l,...,k- 

...,F  .  ,G  . 

n—  l,n  n— l,n 

1. 


1 

}, 


(7.4) 

(7.5) 


From  the  dependency  graph  shown  in  Fig.  7.1,  we  note  that  LRED(i)  and  SRED(i)  will 
always  start  and  complete  at  the  same  time  step  although  SRED(n  —  k  —  1),  the  last 
sequence  for  the  small  block,  will  finish  before  LRED(fc  —  1)  when  k  >  nfa.  In  this  case 
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the  qd  step  of  ^+3,F^+3  . .,  FJ1_l  n  can  proceed  while  the  large  block  is  being 

reduced. 

There  are  at  least  two  storage  schemes  that  can  be  used  for  implementing  the 
bidiagonal  reduction:  cyclic  storage  scheme  and  consecutive  storage  scheme. 

7.3.1  Basic  Procedures 

In  Table  7.1  we  describe  communication  primitives  for  the  node-to-node  commu¬ 
nications. 


Table  7.1.  Communication  Primitives 


Primitives 

Description 

CMMD  Routines 

send(  nodelist;outlist) 

Send  variables  in  outlist  to 
each  node  in  nodelist 

CMMD_send_block 

rec v(i;inlist) 

Receive  variables  in  inlist  from 
Node(i) 

CMMD_receive_block 

s wap  ( i;  inlist;outlist) 

Exchange  variables  in  outlist 
with  inlist  of  NoDE(i) 

CMMD  .swap 

send^and_recv 

( i,j;inlist;outlist ) 

Send  variables  in  outlist  to 
Node (j)  and  receive  vari¬ 
ables  in  inlist  from  Node(z) 
simultaneously 

CMMD.S  end.and  _r  e  cei ve 

Note  that  all  of  CMMD  communication  routines  used  are  blocking  version,  that 
is,  each  node  waits  until  it  finish  sending  or  receiving  the  data  without  proceeding  to  the 
next  executable  code.  This  ensures  that  each  node  carries  out  the  rotation  with  correct 


data  as  we  will  see  in  the  next  section. 
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7.3.2  Cyclic  Storage  Scheme 

Suppose  we  store  the  data  rowwise,  so  that  NODE(i)  contains  all  of  nonzero  entries 

on  the  row  i  such  as  the  diagonal  entry,  and  the  nonzero  entries  created  by  the  chasing 

steps.  Then  we  immediately  observe  that  this  scheme  would  not  give  a  full  parallelism 

among  the  nodes.  For  instance,  at  time  Step  6  in  Fig.  7.1  (sixth  matrix  in  Fig.  6.2), 

G  and  G  cannot  be  executed  in  fully  parallel  fashion  because  Node(2)  has  all 
1,2  3,4 

three  elements  in  the  second  row,  but  (2,3)  entry  is  also  used  in  processing  G„  ,  so  that 
Node(2)  has  to  communicate  with  both  Node(1)  and  Node(3).  In  fact,  any  adjacent 
pair  of  right  rotations  would  cause  similar  difficulties  when  storing  the  data  rowwise. 

However,  this  problem  can  be  completely  avoided  by  storing  the  data  columnwise. 
We  show  this  using  the  following  proposition. 

PROPOSITION  7.1.  Suppose  n  is  even.  Then  if  k  ^  n/2,  the  dependency  graph  for 
Algorithm  6.1  shown  in  Fig.  7.1  can  contain  no  adjacent  pair  of  left  rotations  at  any 
time  step. 

Proof.  We  only  consider  the  reduction  steps  for  large  block  since  the  same  ar¬ 
gument  applies  to  the  small  block.  Let  G.  and  be  the  first  rotations  in 

Lred(z)  and  Lred(i  +  1),  respectively.  Then,  we  see  that  G^  of  LRED(t  +  1) 
is  executed  after  completing  the  rotations  G.  .  .,F.  .  of  LRED(i),  i.e.,  when 

IjIt-  1  I  jIt  i.  1  1 

G.  j  •  °f  LRED(i)  is  executed.  So,  their  indices  differ  by  two.  Since  the  reduction  pat¬ 
terns  proceeds  as  GFFGFG  ■  •  ■  FG,  and  the  indices  for  the  rotations  decrease  by  one 
for  every  pair  of  F  and  G,  the  indices  of  subsequent  rotations  Lred(i)  and  the  ones  in 
Lred(z  +  1)  differ  by  at  least  two.  Therefore,  it  is  impossible  to  have  F |  ^  in  Lred(») 
and  F.  .  in  LRED(f  +  1),  where  \l-j\  =  2  at  the  same  step.  For  Lred(i)  and  Lred(j'), 
where  \i  -  j\  >  2,  the  result  is  more  obvious.  0 
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If  n  is  odd,  we  have  the  same  result  regardless  of  the  value  of  k.  In  fact,  if  n 
is  even  and  k  =  n/2,  it  is  possible  to  have  two  adjacent  left  rotations,  namely,  F, 

•C  1  jK 

and  F ^2  in  the  same  time  step.  But,  the  reduction  step  in  this  case  will  have  the 
following  form: 


However,  since  there  is  no  data  dependency  among  the  nodes,  they  can  be  executed  in 
parallel.  Therefore,  Node(j)  stores  all  of  nonzero  entries  of  j-th  column  including  the 
entries  of  x  and  y,  that  is,  we  partition  the  set  of  nodes  {0, 1, . .  ,,p  —  1}  into 


P  =  {0,1, ...,po—  1},  Ps  =  {p0,p0+  1} 


where  p^q  <  k  <  (pQ  +  1  )q.  Here,  q  =  n/p. 


Node(j'),  j  €  P (,  stores  columns  l  <  k  where 


3  +  1  =  1  (mod  pQ ) 


Node(j),  j  6  P  ,  stores  columns  k  <  l  <  n  where 


j  ~  PQ  +  1  =  (l  ~  k)  (mod  po) 


Example  7.1.  Suppose  p  =  8,  n  =  16,  k  =  11.  Then  q  =  2,  PJ  =  {0,1, 2, 3, 4}  and 
Ps  =  {  5,6,7). 


j 

0 

1 

2 

3 

4 

5 

6 

7 

columns 

1,6,11 

2,7 

3,8 

4,9 

5,10 

12,15 

13,16 

14 

7.3. 2.1  Chasing  Patterns 

Before  we  consider  designing  the  code  for  the  reduction  steps,  we  need  to  catego¬ 
rize  possible  chasing  patterns.  All  of  the  reduction  steps  for  the  large  block  in  Fig.  6.2 
fall  in  one  of  the  chasing  patterns  described  in  Fig.  7.2. 


I  l 


Fig.  7.2.  Chasing  Patterns 
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matrices  and  the  chasing  direction  for  each  block  are  exactly  opposite.  Here,  6’s  represent 
the  diagonal  entries,  a’s  super  or  subdiagonal,  and  d’s  bulges.  Then,  we  see  that  LRED(i) 
and  Sred(i)  have  chasing  patterns,  respectively, 


_ 2U7U _  2(i— 1) 

{/l’Z2’*,3’/4’/3’C-*’Z3’0  ^  . 


Note  that  G ^  2  has  l ^  for  Lred(1)  with  a  undefined  and  l ^  for  LRED(i),  i  >  1.  Similarly, 
G  has  s  for  Sred(1)  with  a  undefined  and  s  for  Sred(i),  i  >  1. 

n— l,n  1  4 

7.3. 2. 2  Host  Program 

The  host  program  distributes  the  data  among  the  nodes,  coordinates  the  order  of 
the  reductions,  and  initiates  the  reduction  process.  The  CMMD  routine  CMMD_distrib_ 
to_nodes  provides  efficient  ways  of  allocating  er.  and  z.  into  the  local  variables  b  and  z 

X  X 

of  Node(z). 

Except  for  the  Lred(1)  and  Sred(1),  the  host  program  initiates  the  subsequent 
reductions  at  every  three  time  steps,  i.e.,  LRED(t)  Sred(w  -  i  +  2)  begin  at  time  step 
3  *  (i  —  1).  This  can  be  done  by  sending  a  signal  to  Node(i)  and  NODE(n  —  i)  as  soon 
as  receiving  the  message  from  Node(*  —  2)  which  just  finished  processing  the  rotation 
F. 

t-2, i-l 

Upon  the  completion  of  Step  1  of  Algorithm  6.1  by  the  nodes,  the  host  program 
calls  the  subroutines  which  will  perform  2x2  updating  or  downdating  and  one  step  of 
the  qd  process  to  complete  the  bidiagonal  reduction.  This  step  of  the  qd  process  is 
completely  serial  unless  k  >  nj 2.  If  it  is  the  case,  Node(£  +  1),  . . .,  and  NoDE(n)  can 
carry  out  the  qd  step  while  the  reduction  on  the  large  block  is  performed. 


7,3. 2. 3  Node  Program 


The  CM-5  allows  a  single  node  program  for  all  the  nodes.  However,  we  need  to 
design  the  code  so  that  a  single  program  can  handle  multiple  chasing  patterns  at  the 
same  time.  Moreover,  each  chasing  pattern  needs  the  data  from  two  or  three  nodes,  and 
these  nodes  must  know  in  advance  what  types  of  operations  to  perform.  For  instance, 
all  patterns  except  for  and  s ^  require  two  nodes  to  do  the  job,  but  and  s ^  need  three 
columns,  and  so  require  three  nodes  to  get  involved  in  the  computation.  To  this  end, 
each  node  keeps  a  local  variable  op  which  is  continually  modified  at  every  time  step. 
The  value  of  op  of  Node (i)  is  determined  by  Node(z  +  1)  for  reducing  the  large  block 
and  by  Node(z  —  1)  for  reducing  the  small  block  as  we  will  see  shortly. 

The  node  programs  for  the  corresponding  patterns  are  given  in  Table  7.2.  A 
unique  value  of  op  is  assigned  to  each  segment  of  the  node  program.  The  value  of 
op  determines  which  operation  each  node  should  perform  at  a  given  time  step,  and  each 
node  executes  only  the  part  of  the  code  marked  by  its  current  value  of  op  . 

The  host  program  ’wakes  up’  the  NoDE(i)  and  NoDE(n  -  i)  always  with  op  = 
1  and  op  =  11,  respectively,  to  start  LRED(i)  and  Sred(z)  because  they  begin  with 
eliminating  x.  and  y  .  .  For  instance,  let  us  consider  Lred(3)  =  {G  ,  F  ,  F  , 

i  n—k— t+1  3,4  3,4  2,3 

G  ,  F  ,G  },  which  has  chasing  patterns  {l  ,1  ,1  ,  l  ,  l  ,1  }.  Then,  Node(3)  receives 
the  value  of  op  =  1,  and  it  immediately  signals  Node(4)  with  op  =  5  to  carry  out  the 
pattern  /  .  Upon  the  completion  of  Z^,  Node(3)  and  Node(4)  increment  their  values  of 
op  by  one  to  continue  on  to  the  next  pattern  Z  .  At  this  point  it  is  not  necessary  for 
Node(3)  to  signal  Node(4)  to  specify  the  types  of  operations.  When  l  is  completed, 
Node(3)  again  increments  its  value  op  by  one  and  signals  Node(2)  with  the  updated 
value  of  op  to  start  the  rotation  F  with  pattern  Z  .  Then,  Node(2)  signals  Node(3) 
with  op  =  7  and  Node(4)  with  op  =  9,  and  all  three  nodes  execute  parts  of  the  code 
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Table  7.2.  Node  Programs  Using  the  Cyclic  Storage  Scheme 


Node(«  +  2) 


Pattern 

Node(i') 

Node(z  +  1) 

/*  op  =  1  */ 

/*  op  =  5  V 

swap(2  +  l;a,  6,x;6,  x) 

swap(i;  6,  x;  a,  6,  x) 

formrot(x,x,c,s) 

formrot(x,x,c,s) 

a  <—  c  *  a 

d  < - 5  *  a 

b  c  *  b 

a  s  *  fc 

d  < - s  *  b 

6  <—  c  *  6 

/*  op  =  2  */ 

/*  op  =  6  */ 

‘2 

formrot(6,  d,  c,  s) 

send(i  +  1 ;  c, 

recv(j;  c,  s) 
applyrot(a,  6,  c,  s,  1) 

/*  op  =  3  */ 

/*  op  =  7  */ 

i3 

recv(z  -f  2;c,  s) 

recv(z  +  2;c,  s) 

6  <—  c  *  6 
d  <—  s*b 

applyrot(6,  a,  c,  s,  1) 

/*  op  =  4  */ 

/*  op  =  8  */ 

swap(z  +  1;  a,  6,  d;  a,  6) 

swap(i;  a,  6;  a,  6,  d) 

i4 

formrot(6,d,  c,  s) 

formrot(fc,  d,  c,s) 

a  <—  c  *  a 

d  5  *  a 

applyrot(a,  6,  c,  5, 1) 

applyrot(a,  6,  c,  s,  1) 

/*  op  =  11  */ 

/*  op  =  15  */ 

swap(i+  l;6,y;  a,6,y) 

swap(i;  a,  6,  y;  6,  y) 

si 

formrot(t/,  t/,  c,  s) 

formrot(y,  y,  c,  s) 

i>  <—  c  *  6 

d  < - s  *  6 

a  <—  s  *b 

b  <—  c*b 

d  <—  s  *  a 

a  4—  c  *  a 

/*  op  =  12  */ 

/*  op  =  16  */ 

S2 

formrot(fr,  d,  c,  s) 

recv(i  +  l;c,s) 
applyrot(a,  6,  c,s,l) 

send(i;  c,  5) 

/*  op  =  13  */ 

/*  op  =  17  */ 

S3 

formrot(a,  d,  c,  s) 

send(i  +  l,i  +  2;c,s) 

recv(i;c,  s) 
applyrot(6,  a,  c,  s,  1) 

/*  op  *  14  */ 

/*  op  =  18  */ 

swap(i  4*  1;  a,  £>;  a,  6,  d) 

swap(i;a,6,d;  a,  6) 

S4 

formrot(6,  d,  c,  s) 

formrot(6,  d,  c,  5) 

applyrot(a,  6,  c,  s,  1) 

applyrot(a,  6,  c,  s,  1) 

d<-s*a 

a  <—  c*  a 

/*  op  =  9  */ 
formrot(a,  d,c,  s) 
send(x,  i  -f  1;  c,  s) 


/*  op  =  19  */ 
rec  v(i;c,$) 
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according  to  the  value  of  op  .  When  this  step  is  finished,  Node(2)  signals  the  host  node 

to  start  Lred(4),  and  it  moves  on  to  the  next  pattern  Z  which  can  be  done  exactly  the 

same  way.  Similarly,  we  process  the  repeated  patterns  Z  and  Z  in  the  same  fashion. 

3  4 

The  node  program  for  the  cyclic  storage  scheme  is  described  in  the  Appendix. 

7.3.3  Consecutive  Storage  Scheme 

In  this  scheme  consecutive  blocks  of  the  bidiagonal  matrix  are  stored  in  each  node. 
We  partition  the  set  of  nodes  exactly  the  same  way  as  in  (7.6)  Using  Matlab  notation, 
Node(j)  stores 


jq+  l:(j  +  1)<7 
(P0-l)?+l:fc 
n:  — 1:  n  —  q  +  1 
n-{p-PQ-  l)?:  —  l:fc  +  1 


if  3  <  PQ  ~  1 
if  j  =  P0  ~  1 
if  PQ  <  3  <  P  ~  1 
if  j  =  P  -  1 


columns  of  B. 

Denote  n  (j)  as  the  number  of  columns  which  the  Node  (j)  has  in  its  memory. 
Then,  we  have 

nc{j)  =q,  i  =  0,...,po-2,po,po  +  l,...,p-2 
nc(P 0“1)  =(!  ~P0)<l  +  k 
nc(P  ~  1)  =  (n-k)-  (p-pQ-  l)q. 

EXAMPLE  7.2.  Suppose  p  =  8,  n  =  16,  k  =  11.  Then  q  =  2,  P  =  {0,1, 2, 3, 4}  and 
P  =  {5,6,7}.  Each  node  has  the  following  columns: 
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j 

0 

1 

2 

3 

4 

5 

6 

7 

columns 

1,2 

3,4 

5,6 

7,8 

9,10,11 

16,15 

14,13 

12 

np) 

2 

2 

2 

2 

3 

2 

2 

1 

Note  that  the  last  n-k  columns  are  stored  in  the  nodes  in  P  in  reverse  order.  This 

5 

ordering  makes  it  possible  for  the  nodes  in  P  and  P^  to  have  identical  node  programs  for 
the  bidiagonal  reductions  on  the  large  and  small  blocks.  From  Fig.  6.2,  we  see  that  the 
bidiagonal  reductions  on  these  two  blocks  are  exactly  opposite,  so  that  the  reduction  on 
the  small  block  can  be  done  by  reversing  the  order  of  the  diagonal  entries  and  y,  reducing 
the  block  exactly  the  same  way  as  the  large  block,  and  again  reversing  the  diagonal  and 
bidiagonal  entries  when  completed.  Therefore,  we  only  describe  the  reduction  steps  (6.8), 
that  is,  only  for  the  nodes  in  P{. 

As  in  cyclic  storage  scheme  NoDE(j)  contains  variables  a[i] ,  b  [i] ,  z  [i] ,  i=l , . . . , 
nc(j)  to  store,  respectively,  the  subdiagonal  or  superdiagonal  entries,  diagonal  entries, 
and  z  .  As  mentioned  before,  the  orders  of  array  elements  are  reversed. 

Following  the  notation  in  (7.4),  Node(O)  initiates  Lred(i),  1,  sequen¬ 

tially.  Then  it  needs  to  communicate  with  Node(1)  to  execute  the  chasing  pattern  / 
of  LRED(g),  reducing  the  entries  of  z.  The  rest  of  LRED(q)  and  starting  LRED(g  +  1), 
which  is  the  responsibility  of  Node(1),  are  done  simultaneously.  We  repeat  this  process 
until  all  of  the  nodes  will  have  finished  their  portion  of  the  bidiagonal  reduction. 

Let  us  call  a  node  which  initiates  a  LRED(i)  for  some  i  at  a  given  time,  the  master 
node.  A  slave  node  is  the  node  which  once  became  the  master  node,  but  now  has  task 
of  chasing  bulges  as  far  as  it  can.  Hence,  at  any  time,  there  can  exist  only  one  master 
node,  say,  Node(j'),  for  some  j  <  p  —  1,  and  j  slave  nodes,  Node(O),  . . .,  Node(j  —  1). 
The  rest  of  the  nodes  stay  idle.  Note  that  Node(po  —  1)  never  becomes  a  slave  node. 

A  pseudo-code  for  the  node  program  is  given  in  Fig.  7.3. 


/*  We  assume  buf  contains  b[q-l],  b[q],  a[q]  ,  z[q].  */ 

/*  Node(O)  is  the  master  */ 
if  myid  ==  0  then 

for  i  =  1  to  numcol-1 
reduce (i) ; 
end; 

send(myid+l,  buf);  recv(myid+l,  buf);  bchase(q,  numcol) ; 
for  i=q  to  k-q  /*  Now  it’s  a  slave  */ 

send(myid+l ,  buf);  recv(myid+l,  buf);  bchase(q,  numcol); 
end; 

/*  Node(idhi)  is  the  master  */ 
else  if  myid  ==  idhi  then 

recv(myid-l,  buf);  redxnode(a,  b,  d,  z,  buf);  send(myid-l,  buf); 
for  i  =  1  to  numcol-1 

reduced);  recv(myid-l ,  buf);  chxnode(a,  b,  d,  buf); 
send(myid-l,  buf); 
end; 

/*  Node(l) ,  ...  ,  Node(idhi-l)  */ 
else 

recv(myid-l,  buf);  redxnode(a,  b,  d,  z,  buf); 
s end (my id- 1,  buf); 

for  i  =  1  to  numcol-1  /*  It’s  now  the  master  */ 

reduced);  recv(myid-l,  buf);  chxnode(a,  b,  d,  buf); 
send(myid-l,  buf); 
end; 

/*  reduce  across  the  node  */ 

send(myid+l,  buf);  recv(myid+l,  buf);  bchase(q,  numcol); 
for  i  *  (myid+l)*numcol  to  k-1  /*  It's  now  a  slave  */ 
send_and_recv(myid-l ,  myid*M,  bufin,  buf  out)  ; 
chxnode(a,  b,  d,  bufin); 

send_and_recv(myid+l ,  myid-1,  buf out,  bufin); 
bchase(q,  numcol); 
end; 
end; 


Fig.  7.3.  Node  Program  for  the  Consecutive  Storage  Scheme 
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Here,  idhi  =  p  —  1  and  numcol  =  n  (j).  The  function  reduce  (i)  performs  the 
following  transformation: 

1  I 


followed  by  bchase(i) 


1  1 


redxnode  carries  out  the  sequence  of  pattern  Z  ,  and  chxnode 
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7.3.4  Computation  Cost 

Let  r  ,r  ,  and  r  be,  respectively,  the  time  required  for  one  multiplication,  one 
MAS 

addition,  and  one  square  root.  Then  the  time  required  for  (6.8)  is 

m  ■=  +  (t  -  l)r.  +  («  -  “)rM  +  E  T,Jk  -  ‘  +  1)  (7.7) 

t=l 

where  T  ,  T  and  T  (i)  are,  respectively,  the  time  required  for  formrot,  applyrot,  and 

g  a  red 

reduce (i).  Since 

T,  =Wtm  +  2ta  <7-8) 

T.  =4r«  +  2rx  P-9) 

TJj)  =(2i-3)r  +(2j-3)r  +(4j-10)rM,  (7.10) 

the  time  required  for  reducing  the  large  block  to  bidiagonal  form  is 

T(k)  =  (16 k2  +  34 k  -  196 )r  +  (4k2  +  8 k-  38 )r  (7.11) 

Moreover,  2x2  updating  or  downdating  requires 


T 

22 


12r  +  r  +  r  for  updating 

MAS  y  b 

9r  +  6r^  +  6r^  for  downdating 


and  one  step  of  the  qd  process  needs 


(7.12) 


T  ,  =  (lOn  —  10k  -  8)r,  +  2(n  -  k  -  l)r 

qd  v  '  M  x  ' A 


(7.13) 
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Therefore,  total  cost  required  for  the  Algorithm  6.1  by  each  node,  is,  on  average, 

v=;<<r(*)+T<"-*)+rM+y>  <7-i4> 

Suppose  k  =  n/2.  Then 

Tb  =1{(2 m  +  T2  +  Td)) 

=  -{(32k2  +  78 k  -  390)r  +  (8fc2  +  18*  -  75)r  +  3 r  } 

p  ivl  A  i i> 

Here  we  took  the  average  for  T  . 

zz 

7.3.5  Communication  Cost 

For  our  analysis,  we  ignore  any  possible  contention  problem  for  simplicity.  Let 

a  and  /?  be,  respectively,  the  startup  time  and  the  data  transfer  rate,  time  per  byte. 
dd 

We  need  at  most  four  real*4  in  single  precision,  so  that  sending  this  message  of  size  16 

takes  a  +4*4/3  seconds.  Then  the  data  in  the  packets  used  in  communication  would 
d  d 

be  of  size  multiple  of  16,  which  makes  the  communication  the  most  efficient  on  CM-5 
[90]. 

7.3. 5.1  Cyclic  Storage  Scheme 

From  Table  7.2,  we  see  that  F  and  G  require  communication  at  every  window 

tj  i  j 

step.  Hence,  from  (7.2)  and  (7.3),  the  communication  cost  becomes 


(7.15) 

(7.16) 


{(5 k  —  9)(ad  +  16/?j)  if  k  >  n/2 

(5 n  —  5k  —  &)(ad  +  16/3p  if  k  <  n/2 


(7.17) 
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seconds,  which  is  impractical  on  a  coarse-grain  distributed-memory  MIMD  multiproces¬ 
sor  due  to  enormous  communication  overhead.  Moreover,  very  few  computations  are 
carried  out  by  each  node  (see  Table  7.2).  Therefore,  a  fine-grain  multiprocessor  with 
fast  communication  ability  would  be  required  to  make  this  storage  scheme  practical  for 
bidiagonal  reduction. 

7.3. 5. 2  Consecutive  Storage  Scheme 

Since  Node(O)  initiates  and  terminates  the  bidiagonal  reduction  for  the  large 
block,  the  total  communication  cost  for  this  scheme  is  determined  by  the  number  of 
send  and  recv  calls  done  by  Node(O).  For  reducing  the  first  q-  1  columns  of  B,  which 
is  done  by  Node(O),  no  communication  is  required.  As  a  slave  node  Node(O)  needs  to 
communicate  k  -  q  +  1  times  with  Node(1)  to  chase  the  bulges. 

Moreover,  when  k  <  n/2,  Node(pq)  makes  extra  effort  to  complete  bidiagonal 
reduction  by  applying  one  step  of  the  qd  process  which  requires  one  more  communication 
step.  When  k  >  n/2,  nodes  in  P*  can  carry  out  the  qd  step  while  the  large  block  gets 
reduced  by  the  nodes  in  P .  Therefore,  the  total  communication  cost  for  this  storage 
scheme  becomes 

(k  -  q  +  2 )(ad  +  16/3rf)  if  k  >  n/2 

{(n  -  k)  -  q  +  1  +  (| PJ  -  l)}(arf  +  16/3rf)  if  k  <  n/2 

seconds.  Here,  IP  I  denotes  the  number  of  nodes  in  P  .  Note  that  T  becomes  worst 
when  k  <  n/2  because  the  reduction  on  the  large  block  is  already  completed  before  the 
qd  process  begins,  and  hence  there  is  no  parallelism  among  the  nodes  for  this  phase  of 
the  algorithm,  requiring  |P  |  —  1  additional  communication  steps. 

S 


T  :=  < 

CON  ' 
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7.4  Timing  Results 

We  implemented  Algorithm  6.1  with  both  storage  schemes  on  CM-5  that  consists 
of  up  to  32  nodes  at  the  Northeast  Parallel  Architectures  Center  (NPAC)  at  Syracuse 
University.  Each  node  of  the  CM-5  is  a  SPARC  chip  which  runs  at  32  MHz  and  delivers 
22  Mips  and  5  Mflops.  There  is  a  64  Kb  instruction  and  data  cache  and  a  16  Mb  memory 
in  each  node  [ill].  In  each  node,  there  are  two  vector  units;  each  vector  unit  is  capable 
of  peak  rate  64  Mflops. 

We  generated  an  n-by-n  random  diagonal  matrix  B  of  the  form  (1.2)  and  a  random 
n- vector  z  with  k  =  n/2,  where  n  =  64,128,256,512,1024.  We  only  implemented  the 
bidiagonal  reduction  part  for  modifying  the  SVD  described  in  Algorithm  6.1.  Computing 
the  SVD  of  a  general  full  matrix  on  a  distributed-memory  multiprocessor  is  described  in 
[71],  and  computing  the  SVD  of  the  bidiagonal  matrix  in  [33]. 

The  execution  time  for  bidiagonal  reduction  using  the  consecutive  storage  scheme 
with  different  matrix  sizes  and  various  set  of  processing  nodes  is  given  in  Table  7.3.  Here, 
we  use  the  following  definition  for  speed-up  achieved  by  our  algorithm 

q_m 

"  T(P) 

where  T(p)  is  time  required  to  execute  the  program  on  p  processors.  Similarly,  we  define 
the  efficiency  of  a  parallel  algorithm  as 

r_  T(l) 

PT(Py 

We  observe  a  linear  speed-up  as  n  increases  with  p  fixed.  When  n  is  small,  we 
cannot  expect  any  speed-up  mainly  because  of  high  communication  cost  compared  to  the 
computation  cost.  For  instance,  when  n  =  64,  more  than  50%  of  the  total  cost  accounts 
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for  the  communication  regardless  of  the  number  of  processing  nodes  (except  for  the  case 
p  =  4).  In  fact,  when  p  =  32,  most  of  time  were  spent  in  communication.  In  this  case, 
since  the  value  of  q  is  only  two,  each  node  has  so  little  to  do,  and  therefore  the  ratio  of 
the  time  spent  in  communication  and  computation  is  quite  large. 

From  (7.18)  since  q  =  n/p,  we  see  that  as  p  increases,  so  does  T  .  However,  in 
general,  q  is  small  compared  to  k.  Hence,  the  difference  in  the  communication  cost  be¬ 
tween  different  p  is  not  significant  as  seen  in  Table  7.3.  In  general,  as  n  increases,  so  does 
q ,  and  therefore,  the  ratio  of  communication  cost  and  total  cost  also  decreases.  However, 
when  p  is  small  and  n  is  large,  significant  part  of  time  were  spent  in  communication 
(p  =  4,  n  =  512, 1024;  p  —  8,  n  =  1024),  and  little  speed-up  was  gained. 

Table  7.4  shows  the  difference  in  execution  time  for  various  k  when  n  =  1024  and 
p  =  32.  Clearly,  we  achieve  an  optimal  speed-up  when  k  is  close  to  n/2.  Note  also  that 
the  bidiagonal  reduction  with  k  >  n/2  is  slightly  faster  than  that  with  k  <  n/2  mainly 
because  when  k  >  n/2,  the  qd  step  can  be  executed  in  parallel  with  the  last  chasing 
step  for  the  large  block. 

We  also  implemented  Algorithms  6.1  and  6.2  using  the  cyclic  storage  scheme. 
However,  it  was  embarrassingly  slow,  and  no  speed-up  was  gained  in  any  case.  As  we 
have  shown  in  the  previous  section,  this  storage  scheme  becomes  totally  impractical  when 
n  >  p  due  to  high  communication  overhead  caused  by  severe  contention  problems. 
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Table  7.3.  CPU  Time  (in  sec)  for  the  Bidiagonal  Reduction 


Table  7.4.  CPU  Time  (in  sec)  with  Various  k  (n  =  1024,  p  =  32) 
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Chapter  8 


Conclusion 


We  have  presented  efficient  algorithms  for  modifying  the  TSODs,  and  shown  that 

2 

the  TSOD  can  be  updated  and  downdated  in  0(n  )  flops  in  a  manner  that  preserves  its 
structure.  The  backward  error  analysis  and  perturbation  theory  show  that  the  proce¬ 
dures  satisfy  a  blockwise  stability  property.  Thus  if  our  interest  is  in  separating  the  two 
subspaces  associated  with  the  large  and  small  singular  values,  we  will  obtain  answers 
that  are  as  good  as  can  be  expected.  The  use  of  this  perturbation  theory  shows  that 
we  achieve  more  accuracy  in  the  singular  values  and  more  orthogonality  in  the  singular 
vectors  that  result  from  our  update  procedures.  Our  numerical  tests  show  some  im¬ 
provement  in  the  accuracy  of  downdated  singular  values  using  our  algorithm  instead  of 
general  chasing. 

Our  approach  to  modifying  the  ULVD  is  particularly  promising.  It  is  simple  to 
implement  for  both  updating  and  downdating,  and  preserves  the  rank-revealing  struc¬ 
ture  often  without  the  deflation  process  to  compute  the  numerical  rank  of  the  matrix. 
Moreover,  one  can  efficiently  track  the  size  of  each  block  of  lower  triangular  part  of  factor¬ 
ization  for  an  accurate  monitoring  of  the  condition  of  downdating  problem.  Furthermore, 
data  independence  among  the  blocks  makes  the  algorithms  parallelizable. 

We  also  have  given  algorithms  for  rank-one  updates  and  downdates  of  the  SVD 
and  partially  reduced  bidiagonal  forms.  It  has  also  been  shown  that  these  algorithms 
satisfy  a  blockwise  stability  criterion  that  has  not  been  shown  for  other  algorithms.  The 
algorithms  proposed  to  allow  the  user  to  specify  a  tolerance  between  large  and  small 
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singular  values,  and  the  separation  between  the  associated  subspaces  is  preserved  by  the 
algorithm. 

Finally,  we  presented  an  efficient  method  for  modifying  the  SVD  in  parallel  to 
establish  its  practical  value  in  real  time  applications.  The  algorithm  preserves  the  block 
structures,  maintaining  the  efficiency  of  parallel  procedures  as  well.  The  consecutive 
storage  scheme  outperforms  the  cyclic  storage  scheme  in  bidiagonal  reduction  due  to 
high  communication  cost  of  the  latter  scheme.  The  experiments  show  the  efficiency  in 
using  the  processors  was  slightly  over  60%  for  the  consecutive  storage  scheme. 

Although  the  entire  thesis  is  devoted  to  the  problem  of  modifying  the  TSOD,  there 
are  a  number  of  unresolved  issues  and  problems.  We  suggest  a  few  in  the  following: 

•  The  stability  of  the  CSNE  approach  taken  when  downdating  is  not  possible  is  de¬ 
termined  only  by  experimental  results.  In  particular,  one  needs  to  know  how  good 
x  in  (6.17)  really  is.  A  rigorous  error  analysis  should  be  performed  by  extending 
such  results  as  those  due  to  Bjorck  [18]. 

•  Parallelizing  the  ULVD  procedures  is  much  more  challenging  than,  the  SVD  pro¬ 
cedures.  Although  the  reduction  on  the  large  and  small  blocks  can  be  executed 
simultaneously,  it  is  difficult  to  enhance  parallelism  within  the  block.  As  suggested 
in  [103]  for  parallelizing  the  URVD,  we  may  require  a  fine-grain  MIMD  architecture 
for  efficient  implementation  of  systolic  arrays. 

•  A  procedure  described  by  Gu  and  Eisenstat  [56]  is  a  promising  approach  to  modify 
the  singular  vector  matrix  after  modifying  the  SVD.  They  use  an  adaptive  version 
of  one-dimensional  fast  multipole  method  [29].  However,  an  efficient  parallel  imple¬ 
mentation  for  this  acceleration  method  is  not,  as  yet,  available  although  a  parallel 
procedure  for  non-adaptive  version  has  been  developed  [45].  Together  with  our 
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parallel  bidiagonal  reduction  procedure  described  in  Chapter  7,  and  parallel  al¬ 
gorithms  for  computing  the  SVD  of  bidiagonal  matrices  [33,  65,  107],  the  SVD 
algorithm  described  in  Chapter  6  can  be  implemented  fully  in  parallel. 

•  It  would  also  be  interesting  to  extend  algorithms  for  modifying  the  ULVD  to  those 
for  modifying  the  ULLVD  for  two  matrices  [75]  as  an  approximation  to  the  gen¬ 
eralized  SVD  [84,  120].  We  also  need  to  analyze  the  stability  and  complexity  of 
the  algorithms  for  modifying  the  ULLVD  and  derive  the  error  bounds  on  the  sub¬ 
spaces  computed  by  the  ULLVD  compared  with  those  by  the  generalized  SVD.  The 
Estimator-Correlator  array  processor  [108,  109]  that  can  implement  the  estimator 
kernel  using  the  ULVD  and  the  inverse  noise  kernel  using  the  ULLVD  would  be  a 
practical  application. 
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Appendix 

Node  Program  for  Bidiagonal  Reduction 
with  Cyclic  Storage  Scheme 


c  This  node  program  will  reduce  the  arrow-head  matrix  to 
c  the  bidiagonal  matrix  using  cyclic  storage  scheme 

c 

c  Variables  used: 

c  abdz  —  contains  values  of  a,  b,  d,  and  z  in  that  order 

c  myid  —  node  id 

c  idpl  —  myid  +  1 

c  idp2  —  myid  +  2 

c  idml  —  myid  -  1 

c  idm2  —  myid  -  2 

c  hid  —  host  id 

c  cs  —  angles  for  the  rotations 

c  inmgl  —  processing  the  large  block 

c  inrng2  —  processing  the  small  block 

c 

special  case  —  first  two  rotations  in  each  block 

100  continue 

again  =  .true. 

50  continue 

ret  =  CMMD_receive_block(src,  tag,  op,  isize) 

150  continue 

go  to  (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19),  op 

1  ...  op  *  1 

2  ...  op  =  2 

3  continue 

ret  =  CMMD_send_block(idpl,  tag,  7,  isize) 
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ret  =  CMMD_send_block(idp2,  tag,  9,  isize) 

ret  =  CMMD_receive_block(idp2,  tag,  cs,  angsiz) 

abdz(3)  =  cs(2)*abdz(2) 

abdz(2)  =  cs(l)*abdz(2) 

op  =  4 

it  =  it  +  1 

go  to  200 

4  continue 

ret  =  CMMD_send_block(idpl,  tag,  8,  isize) 
ret  =  CMMD_swap(idpl ,  abdzt,  bufsiz,  abdz,  bufsiz) 
call  formrot(abdzt(2) ,  abdz(3) ,  cs(l),  cs(2)) 
abdz(l)  =  cs(l)*abdz(l) 

call  applyrot (abdzt ( 1) ,  abdz(2) ,  cs(l),  cs(2),  1) 
op  =  3 
it  =  it  +  1 
go  to  200 

5  ...  op  =  5 

6  ...  op  =  6 

7  continue 

ret  =  CMMD_receive_block(idpl ,  tag,  cs,  angsiz) 
call  applyrot (abdz (2) ,  abdz(l) ,  cs(l),  cs(2),  1) 
it  =  it  +  1 
go  to  100 

8  continue 

if  (inrngl  .and.  it  .eq.  maxit)  then 

ret  =  CMMD_send_block(hid,  tag,  20,  isize) 
end  if 

ret  =  CMMD_swap(idml ,  abdzt,  bufsiz,  abdz,  bufsiz) 
call  f ormrot(abdz(2) ,  abdzt(3),  cs(l),  cs(2)) 
abdz(3)  =  cs(2)*abdzt(l) 

call  applyrot (abdz ( 1) ,  abdzt(2),  cs(l),  cs(2),  1) 
it  *  it  +  1 
go  to  100 

9  continue 

call  formrot(abdz(l) ,  abdz(3),  cs(l),  cs(2)) 
ret  =  CMMD_send_block(idm2,  tag,  cs,  angsiz) 
ret  =  CMMD_send_block(idml ,  tag,  cs,  angsiz) 
it  =  it  +  1 
go  to  100 
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10  ...  terminate 

go  to  999 

200  continue 

if  (again)  then 

again  =  .not.  again 
go  to  150 
else 

if  (myid  .ne.  0)  then 

ret  =  CMMD_send_block(idml,  tag,  op,  isize) 
else 

ret  =  CMMD_send_block(hid,  tag,  op,  isize) 
end  if 
go  to  100 
end  if 

999  continue 

.  similarly  for  the  small  block 


